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VORTEX METHODS FOR SEPARATED FLOWS 

By PHILIPPE R. SPALART 1 


Summary 


1 h, numerical -olulum of the Euler or Navier-Stokes equations by Lagrangian -vortex 

‘ : S “ I— M-und is presented in a„ elemenllry £h 

and eludes the relationship with traditional point-vortex studies, the convergence to 

, V !';' ll " J C,| ' JaUo " s ’ and ,he essential differences between two- ,„d 

i'we ;l3: 1,1 cxien,iin ’ 1 ,w mcihod - - — ««* 

einnl'e "'!' •"** «*■» "*><*>' » kept to a minimum and more 

3ies ’’me "" A “i, ar ™ 0f ” PCT,iK ' namd >' t wo *dimensional flow, around 

I In I bodies. When soud walls are present, complete mathematical results are not available 

aim one must adopt a more heuristic attitude. The imposition of inviscid and viscous 

wXmT-iT 7T ' V, |! ™ hmU> "“«’*"»• «■-* vortices and the creation of 

i< tv along solid wall, are examined in detail. Methods for boundary-layer treatment 
* 1 *' question of the Kutta condition are discussed 

Practical a,, lects and tips helpful in creating a method that really works are explained. 

, , !"" of »he method and the assessment of accuracy, vortex- 

<• « pio lies, tiiiie-inardnng schemes, numerical dissipation, and efficient programming 
Operation counts for unbounded and periodic flows are given, and two algorithm designed 
"P the calculations are described. 8 

Calcidations of flows past streamlined or bluff bodies are used as examples when appro- 

i f i ' ! -7 ,'t f "T T" g lams ' ,h ' 5,ar,i11 * a '"t dynamic stall I,, 

at , - d Hit,! 0 A ” a - W< "' '"‘"' si ° nal cas, ' a<le ' a 'nnlti-elenient airfoil, and an attempt 
at pr. dieting the drag crisis of a circular cylinder. H 
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1. Introduction 

Vortex methods in general were thoroughly reviewed by Leonard (1980, 1985), and it 
would be difficult to improve on these articles. In the present notes we intend to cover the 
basics of vortex methods and then to discuss in more depth two subjects: the interaction 
with solid walls, and the practical aspects of programming and using vortex methods. 
These subjects were of less interest to Leonard, but are crucial in engineering applications. 
I herefore we attempt to present, useful and sometimes original material in these areas. 
We shall also mention some recent contributions, published after Leonard’s articles. The 
sections which contain advanced material and could be omitted for a first reading are 
indicated by stars*. 


1.1. Example: flow past a multi-element airfoil 

V\e propose to start by showing a calculation of a separated flow via the vortex method to 
illustrate its main features and, we hope, render it attractive to the reader. We shall often 
describe the method by comparison with grid-based (finite-difference or finite-element) 
methods, since these are more widespread. 

In 1986 the Boeing Commercial Airplane Company contacted the author, asking for 
computations of the flow past an airfoil in landing configuration, with a leading-edge slat 
and a douhle-slot ted trading-edge flap. Boeing was especially interested in the variation 
of the lift with Reynolds number. The configuration was a severe test for numerical meth- 
ods due to the complex multiply-connected shape and the importance of viscous effects, 
including separation. 

Results were obtained in less than a week. In spite of the many concave and convex 
sharp corners present, no smoothing or other alteration of the shape was necessary. This 
illustrates the first advantage of the vortex method: it is grid-free. With a finite-element or 
especially a finitodiflerence method the grid generation for such a shape would he difficult 
and time consuming, and one may have to settle for rather distorted grids which degrade 
the accuracy. With the vortex method the only precaution taken was to make the spacing 
between vortices and the time step small enough to accomodate the narrow gaps between 
elements 2 and It, and 3 and 4 of the airfoil. 

figure 1 shows the airfoil, the vortices (more accurately, the centers of the vortex blobs), 
streamlines, and the individual (hollow arrowheads) and overall (filled arrowhead) force 
vectors at a given time. Observe the recirculating bubble behind the slat, the acceleration 
of the fluid through the slots, the separation at sharp corners and on the top surface 
of element 1, and the irregular motion in the w'ake. The computed forces were in good 
agreement with experiment. 'The figure illustrates the Lagrangian character of the method: 
the solution procedure consists in tracking vortices which move with the fluid, rather than 
updating quantities at fixed grid points. It also illustrates the second advantage of the 
method: the vortices carry all the information ami are needed only in a narrow region 
near the body and in the wake. In fact the calculation used just 1300 vortices (roughly 
the equi valent of a single 36 - 36 grid). The streamlines were computed only for display 
purposes; for the purpose of advancing the equations in time, computing the forces etc., 
one only needs to compute the velocity of the vortices themselves. 

Additional advantages of the method are first that the boundaries of the computed 


Figure 1. Flow past a multi-element airfoil. ••• vortices; 1 forces; — streamlines 


domain are at infinity which removes all of the problems associated with computations in a 
truncated domain (the effects of the boundary conditions on accuracy and on stability) and 
second that there is low numerical dispersion. The method can transport flow structures 
(groups of vortices) at any velocity without deforming them or dissipating them as a grid- 
based method may do There is no CFL number. The question of time accuracy will be 
addressed in much more detail later. 


Among the disadvantages are the fact that the flow had to he treated as incompressible, 
although fairly high local Mach numbers may he reached in the slots even at the landing 
speed. Th»* calculation also required the storage of a 249 x 249 full matrix and about 
. r > x 10 floating-point operations to establish the flow, starting from rest, and generate a 
time sample long enough for averaging. Both figures are much larger than what a 36 x 36 
grid calculation would require. Thus, a vortex carries much more information than a grid 
point, hut it is also much more expensive to maintain. This is because all die vortices 
interact, so that the operation count is of order A' 2 , where i\ is the number of vortices. 
The competing methods often have operation counts of order N or AMog(A'). 

to answer Boeing’s question, a Reynolds number dependence was indeed predicted. As 
the choid Reynolds number was changed from 2.3 x 10* to 1.2 x 10 7 with all the other 


parameters (physical and numerical) exactly the same, the lift coefficient at 8° incidence 
in* teased by about Ibis was caused by a difference in the boundary-layer behavior 

on the upper surface of element 4. At the lower Reynolds number, the boundary layer 
always separated, whereas at the higher Reynolds number it could transition and remain 
attached part of the time, thus increasing the circulation around element 4 and therefore 
the overall lift. This Reynohls-nuinber trend was consistent with the common wisdom. 
Irouiially. Boeing subsequently informed the author that experiments had predicted an 
opposite trend for this configuration, aU that the situation was rather confused. This 


study demonstrates the power of a well designed vortex method in dealing with complex 
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geometries, and how challenging and unpredictable Reynolds-number effects can be to 
numericians, experimentalists and designers. 


1.2. Governing equations 

Complete discussions of these equations can be found in textbooks, e.g., Batchelor 
(1907). The fundamental equations arc the incompressible Navier-Stokes equations, 


V.U = 0 

V. 4-U.VU -Vp+ t,V 2 u. 


(la) 

(lb) 


Here U is the velocity vector, V is the gradient operator, t is the time, p is the kinematic 
pressure, and /' is the kinematic viscosity. Bold-face letters like U denote vector quantities, 
and the symbol = indicates a definition. Equation (la) is the continuity equation and (lb) 
the momentum equation. One obtains the Euler equations by setting v — 0. The vorticity 
u> is the curl of the velocity, 


w = V < U. 


( 2 ) 


In an unbounded domain with the fluid at rest in the far field, one has the “Biot-Savart 
law”, which is the inverse of (2) by a Green ’s-function approach (when (la) holds) and 
explicitly gives the velocity in terms of the vorticity: 


U(x) 


i r u>(x') x (x - x') 
4n J |x - x'j 8 


dx'. 


( 3 ) 


This equation only expresses the kinematics of the flow. The integral, taken over the whole 
two- or three-dimensional domain, is responsible for the nonlocal interactions between 
vortices. By taking the curl of (lb) one gets the vorticity equation, 


4 U.Vu) ui.VU 4 i'V 2 u>. (4) 

Notice that the pressure term is eliminated, but that a new term, u.VU, appears which is 
responsible for vorticity stretching and rotation. In this term, VU can be replaced by its 
symmetric part (VU | VU 7 )/2 (the deformation tensor), which may be advantageous in 
some numerical methods (Oantaloube k Huberson 1984). 

I he circulation of a velocity field U along a closed contour is defined by the line integral 

r -• f v.di ( 5 ) 

Tins is an extremely useful concept because Kelvin’s theorem (Batchelor 1967, p. 273) 
shows that in the absence of viscosity, the circulation around a contour that follows the 
fluid is conserved: 1)Y / Df 0. 

In two dimensions the equations are much simpler. If the flow is in the (x,y) plane the 
vorticity only has a component in the ; direction, w. =• dv/dr - du/dy (which we’ll denote 
simply by oj). The Biot-Savart law becomes 


U(x) 




uj(x') dx' 


( 6 ) 



Invariants of the motion 

.MlloT' S °!"7 f f ) , ' OU# 1 d ™" m * ,,umeric<J ■"‘‘Hod is its conservation of at least 

f( . appropna e tn . egpil quantities. In the present context, the most useful form 

• . . ^, 9G ’ p> jl ‘), Wu k Sankar (1980) or Ting (1983); the arguments relv 

on m igrations by parts but are subtle, because the velocity decays slowly (Algebraically^ 

situatfoiA ’inThreeT" W ^ W VOrtldty decays fast ("ponentially) which is the usual 
-ituation. In three dimensions the important, integrals are 


u> dx, 


x x u) dx, 


J X x ( x Xu;) <7x, 


u>(x).u-(x') 


dx dx'. 


energy resnAcUv \° w'^ 1,ucar "^omentum, angular momentum, and kinetic 

these YiileuK Is I i I e “““'‘"I ‘ ilal llle ™ rtid, > drra y s <»*l enough in the far field for 
s , v r q T N< "'' " ,C “*! rartid ‘.V (9) is trivial unless * decays 

1,IV laSU ' r ’ '•»» r . one can easily show that this integral 

n.i .v zuo use .ie - 0). In two dimensions the corresponding quantities are 
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Here, the total vorticity (13) has no reason to be zero. 

In unbounded inviscid flows, all these quantities are conserved. In unbounded viscous 
flows, only total vorticity and linear momentum (9, 10, 13, 14) are conserved. Formulas 
giving the rate of viscous decay of angular momentum and kinetic energy can be found 
(e.g., Ting 1983). When there is a solid body in the fluid it can alter any of the integrals 
(9-16). However, the variation of (9) and (13) (with the integrals restricted to the fluid 
region) is known once the motion of the body is (Wu & Sankar 1980). 

The invariants are often invoked when constructing or evaluating methods, but the 
actual importance cf exactly conserving quantities such as (9-16) in a numerical method 
is a matter of debate . The invariants are integrals over the whole domain, and contain i o 
information about the details of the solution. Many successful methods violate some of 
the invariants, especially once time-integration errors are included, and one can construct 
methods that conserve the invariants but are not very good or even are inconsistent. The 
common argument t hat energy conservation prevents the calculation from “blowing up” is 
rather technical, and in any case does not directly apply to vortex methods. 

The following intuitive argument seems more relevant. The conservation of an invariant 
means that over the whole domain the local errors exactly cancel each other in some 
measure (e.g., |x 2 |c/x for (15)). From a statistical point of view, it is then reasonable to 
expect that, over an intermediate domain, large enough to contain many vortices or grid 
points, the local errors tend to cancel each other instead of accumulating. One then expects 
much better accuracy on intermediate scales than on small scales. This is probably true 
in many vortex calculations, especially complex ones such as in Fig. 1. One can tolerate 
significant errors in the irregular small-scale motion, as long as they do not affect the 
quantities of interest such as the forces and moments on the bodies. 

1.1). Constraints on vector fields* 

Vortex methods and related methods focus on derivatives (or combinations of deriva- 
ti'-es) of the primary variables. Using a collection of computational elements (e.g., vortices) 
one can represent an arbitrary vector field. This can become i liability, because many of 
the vector fields involved are not arbitrary; they satisfy constraints. For instance, if the 
vector A is the gradient of a scalar, say A Vd>, then V '< A = 0. Now if one represents 
the components of A separately with two collections of computational elements there will 
in general be an inconsistency because dA r /dy may not be equal to dA y /dx, and if they do 
difier there is no scalar <f> such that d0/d.r = A x and defr/dy = A y . For instance Anderson 
(1985a) proposed a Lagrangian method for the transport of a gradient and avoided this 
question by solving a Poisson equation for <p: V 2 (b = dA x /dx + 3A y /dy , instead of the 
original equations. For the vorticity, the constraint, is that it should be divergence-free: 
V .w 0. As we shall see this also caused problems, and again some methods had recourse 

to a Poisson equation (V 2 U — — V x u>). 

A more direct response to the question of the constraints is to “build them” into the 
discretization. For instance the Contour Dynamics method (Zabusky, Hughes & Roberts 
1979) can be regarded as based on computational elements that “carry” Vu>, the gradient 
of the 21) (scalar) vorticity, but the elements are automatically arranged on chains so that 
the vector field that purports to be Vu: is indeed curl-free. The equivalent for the vorticity 
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is the filament method, in three dimensions. Such “chained” methods are inconvenient in 
terms of bookkeeping (note how the Contour Dynamics method was never applied to more 
than a few chains), but are more attractive intuitively. This will be discussed again in the 
context, of the 3D methods. 
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2. Point-vortex methods 


2.1. Basics 

Point- vortex studies have a long history. Interest in them is justified by the fact that 
die dynamics of point vortices provide exact albeit singular solutions of the incompressible 
equations in two dimensions. Assuming one can perform the time integration either 

n T! ■ y , with hig !; accuracy ’ ° ne can generate n ° ntriviai Euier -j^ions 

with a small computational cost and very little memory 

v „?"f d ' r * ° !N "*«“* P“i«ons 1,JV, with circulations T,. The 

vort icily distribution is 3 

N 

u,(x) = £ r , 6 ( x~ Xj ) (17 ) 

where S is the two-dimensional Dirac distribution (6 is infinite at the origin, 0 elsewhere 
(6) we get fsu," 1 ' ^ ^ ^ & C ° 11CCti ° n ° f “ SpikeS ” ° f V ° rtirity ' Insertin g ( 17 ) int ° 

= as) 


*= i 


t verify that the vdocity fie,d w * ^^4^ is 

only the l^'vOTfex 6 ^ P ° mtS **’ ^ ^ * C,rcuIation Fk around an y contour enclosing 

The other ingredient which expresses the dynamics, is Kelvin’s theorem. If we make 
ach vortex follow the fluid and imagine a small contour around it we get tie result that 
the circulation I j of each vortex is conserved. Therefore: 


dr, 

di 

dx 

di 


~ 0, 


r = U(xj,/). 


(19) 

(20) 


T he formulation is now complete. The unknowns are Xj and T,, j = \ N. The time 
cvo Ution IS given by (19) and (20), with (18) giving U in (20). Note that (18) is singular 
when =-- x*.; tins contribution is simply omitted from the sum. In other words there is 
no se -induced velocity for the vortex. The circulations are constant, and the only task is 
to move the vortices with the fluid. 

An early example of point-vortex studies is Rosenhead’s (1931). As an exercise one can 
piogram the equations (using one’s favorite explicit time-integration scheme for now) for 
wo vor ices with various circulations of the same sign or of opposite signs. The vortices 
should orbit around their centroid (T,x, 4 T 2 x 2 )/(r ] 4 T 2 ) unless their circulations are 
exactly opposite, in which case the pair moves in a straight line. 

Lven though (1 1 ) looks like a simple linear superposition of solutions, it is not, because 
the vertices interact through (18). This is how the nonlinearity of the Euler equations 
elects he vortex equations. Note also that the equations (18,20) can be derived from the 
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Hamiltonian H _ ^ I jT* log | Xj x*. |/ 47r with I'jXj and yj as the canonically conjugate 

coordinates and momenta (Aref 1983). H corresponds to (16) with the self-energy of each 
vortex (which would be infinite) omitted; in that sense it is the “interaction energy”. 
Naturally, H is conserved by the dynamical equations. 

To bo rigorous, the point-vortex method can only represent very special vorticity distri- 
butions, namely superpositions of Dirac distributions as in (17). These Dirac distributions 
are singular, and so is the velocity field they induce: it tends to oo like l/|x - xt| as x 
approaches x k . Such flows are not realistic in the sense that they have unbounded veloci- 
ties, infinite kinetic energy, et c. However, if one considers flows with small but nonsingular 
vortices the point-vortex method is still useful as we show below. 

An example of a realistic thin vortex is the “spreading line vortex” which is a classical 
exact solution of the 2D Navier-Stokes equations (Batchelor 1967, p. 204). The vorticity 
is given by 


w (x, t ) - 


r 




( 21 ) 


4nut r ' 4 ui 

the vortex being centered at the origin (x, ==. 0). As the product ut tends to 0, (21) 
approaches the Dirac distribution (17). j’he velocity induced by this vorticity is 



As ut tends to 0 this velocity tends, point by point, f.. the one given by (18). Thus if 
the distance between vortices is large compared with their radius (v/nd) the equations 
(IS, I!), 20) will not be affected. More generally, if we have small patches of vorticity 
(not necessarily of the type of (21)) of extent c much smaller than the distance d between 
patches the interaction of the vortices through (IN) will not be affected. The patches also 
have internal dynamics; they deform, and spread if viscosity is present. However it can 
be shown that the centroid of the vorticity in the patch is not affected by the internal 
dynamics (convective or viscous) (see (14)). Therefore if x y represents the centroid of the 
|>h tcli, llir sdf-imluced velocity is still exactly zero. 

The conclusion is that as the ratio t/d tends to zero the motion of the patches tends to 
that of point vortices. This was shown rigorously by Marchioro & Pulvirenti (1983) for 
shoit times, and it is reasonable to expect that it holds for all times, except maybe for 
some very special initial conditions. For the gravitational law and well-separated planets 
the analog is the familiar point-mass approximation. Saffman & Meiron (1986) present a 
different argument, showing that point vortices provide a weak solution of the vorticity 
equation (7) with u ~ 0 (although according to C. Greengard & Thomann (1988) this is 
controversial, and depends on one’s interpretation of how the Dirac distribution acts on 
a singular function such as (18)). Thus the point-vortex method is not strictly limited 
to vorticity fields of the type of (17); it is also a good approximation for the motion of 
small, well-separated patches of vorticity. This is still too restrictive for most practical 
applications. Consequently the recent work with point vortices tends to focus on math- 
ematical issues like the Hamiltonian character of the equations, the existence of steady 
< on ligui at ions, integrahle and chaotic motion (Aref 1983), the statistical properties, and 
mixing of the fluid surrounding the vortices. A more elaborate model for well-separated 
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patches is the “moment model” of Melander, Zabusky & Styczek (1986), which includes 
an approximation of the inviscid internal dynamics. 

2.2. Difficulties with vortex sheets * 

Given sufficient computing power it is natural to use more and more vortices to obtain a 
finer and finer description of a given vorticity distribution, just like one refines a grid. This 
idea was first applied to vortex sheets. Such flows are not as strongly singular as point 
vortices: the vorticity is still singular but only like a one-dimensional Dirac distribution, 
the velocity is discontinuous but finite. In fact Rosenhead’s work was aimed at a vortex 
sheet, but he did the computations by hand and could not use a large number of vortices. 
When this became possible, the outcome was a disappointment: smooth sheets could not 
be maintained, and the vortices rapidly started moving chaotically (BirkholT 1962). The 
time integration was accurate, and the situation did not improve as more vortices were 
used (Moore 1979). Ostensibly, the method was not converging. 

The consensus now is that the wrinkling of the sheets discretized with point vortices was 
not just a numerical phenomenon, because the vortex-sheet problem is not well posed: in 
general the exact solution does not remain smooth even if the initial condition is. This 
is so because waves of all scales are unstable, and the shortest waves grow the fastest. 
Krasny (1986) and other researchers have shown that a periodic vortex sheet disturbed 
by a single sine wave remains smooth for a finite tune, and then develops a singularity 
(infinite curvature). Krasny explains that if the initial condition is analytic the amplitude 
of the short waves is exponentially small, so that even though they grow very fast they 
do not contaminate the whole solution until some finite critical time t c . Krasny showed 
that the point-vortex method gave a good solution up to time / r , but only if round-off 
errors were reduced to a low enough level. He suggested using higher-precision arithmetic, 
and also proposed a filtering technique which is very powerful, but should be used with 
exquisite care. 


2.3. Application to smooth Euler solutions* 

1 he point-vortex method is not considered a viable numerical method to converge to 
smooth solutions of the 2D Euler equations. However, the reason is not that it would fail 
in a blatant way, as it seemed to do with vortex sheets. The main reason the point- vortex 
method ;s out of favor is that convergence proofs have been given for vortex-blob methods 
(!|3) an d demand that the blobs overlap, and overlap more and more as they become 
denser. I his is incompatible with point vortices. There is the question of what measure 
of the error one is using. If one takes the velocity or vorticity fields, there is no way the 
point-vortex solution could converge in any of the familiar norms, e.g,, since it is 
even bounded. On t he other hand if one only considers the trajectory of the vortices aim 
their initial arrangement is regular, they behave well at least for some time, which implies 
that the global measures like the moments of the vorticity and the Hamiltonian (13-16) 
behave well too. It may be that the trajectories (and even the velocity field if one takes 
an exotic norm like a Sobolev norm with a negative index) converge to the exact ones 
at least up to some critical time. This is just a conjecture. 
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3. Vortex-blob methods 


3.1. Motivation and basics 

Initially a strong motivation to switch from point vortices to nonsingular (or less- 
singular) elements than was the failure of the method on vortex sheets (Chorin & Bernard 
1973, Kuwahara & Takaini 1973). Now, the emphasis is on smooth Euler solutions, for 
which the blob methods have been shown to converge, even with high orders of accu- 
racy (Hald 1979, Beale fc Majda 1982). We present the method in two dimensions. In a 
vortex-blob method the Dirac distributions b in (17) are replaced by bounded functions 

b n : 

N 

^( x ) =- ^ ^ 6 a (x - Xj). (23) 

j. i 

i he b a function has ai: integral of 1, and the length scale <r (the “core radius’ 7 ) is a 
measure of the spread oj b 0 . In almost all cases, b n is radially symmetric and its shape is 
independent of <r i.c\, 

c / . 1 r/ lx - X,*|v 

^(x - Xj) = —■ /(- — -— ) (24) 

for some nondimmsional function /, usually bell-shaped, and satisfying 

2?r f rf(r)dr 1. (25) 

./(» 

One now has a “realistic” velocity field, hounded and with locally finite kinetic energy. 
The function / is small away from the origin, and as a - (), the blobs resemble point 
vortices. The Gaussian, as in (21 ), is a possible choice but may not be the most judicious 
in terms of computational speed. Equation (IS) becomes 


U(x) 



,v 

Vr 

h i 


X 

|X 


Si , /. ( x x ‘j) 

Xjtr rr 


(26) 


where F is defined by 

F(r) _ 2tt f rf(r')dr f . (27) 

Jo 

Figure 2 gives an example of the velocity given by (26, 27), compared with the point- vortex 
velocity ( 1 ^)- 1 hey agree except for |;r|/<r of the order of 1 or less, and the blob-induced 
velocity is smooth at the origin. Finally the stream function is given by 


V’(x) 


iVr,/-r 


!x 


X* 


2tt 


h i 


(2S) 


with F defined by 


r(r') J , 
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Figure 2. Azimuthal velocity induced by vortices. — vortex blob; - - - point 
vortex 


Equation (19) is unchanged; circulation is conserved. In most methods (20) is also 
unchanged, i.e., the blob moves with the velocity at its center x } . Leonard (1980) discusses 
the use of the velocity averaged over the blob (weighted by 6„) instead of the velocity at 
the center. In that case, (27) does not hold and the relation between / and F is slightly 
more complex; the alteration of the velocity is 0{(r 2 ). Averaging the velocity is sensible 
since the most meaningful interpretation of x is as the centroid of a small vortex patch 
(§2./). Leonard also shows that the kinetic energy of the flow (16) is conserved only with 
the weighted scheme. On the other hand this energy strongly depends on the function 
/, which has little physical meaning, and given a function F one can always define a 
Hamiltonian which approximates the kinetic energy, and will be conserved. The weighted 
scheme has not received much use. 

One could consider using a different core radius n } for each vortex; however this is ill- 
advised because it prevents the method from conserving even the most basic invariants, e.g., 
the linear momentum (14) (Leonard 1980). One can of course symmetrize the interaction 
by taking an average of (Tj and cr*., which restores the conservation of (14), but this is 
arbitrary and seems inconsistent with the concept that the exact value of rij is important. A 
much better response to the objection that (14) is not conserved would be to use Leonard’s 
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.' C,ght ° d SC,,C, ' ir - S,,n ’ tho consc,lsus ls t>'at ilii- calculations arc, and should be, quite 
insensitive to the function / and the value of it. The basic method defined by (19, 20, 26) 
with uniform it (independent of j or /) is by far the most used. 

■i.'i. Convergence in theory 

The accuracy and convergence have been discussed by Hald (1970), Leonard (1980), and 

en e , ajda (1982), among others. Leonard focused on the consistency of the method 
w u reas Hald and Beale it. Majda obtained complete convergence proofs including stability. 
Both approaches produced “moment conditions” which indicate the order of accuracy of 
I ne method according to whether some integrals of the function / equal 0. 

Beale k Majda split the error in the velocity field at some time / into two components, 
one due to the inexact, position of the vortex centers (compared with the motion of the 
points in the exact solution), and the other due to the discretization of the velocity by 
M,). I ne two components naturally interact, as / increases. Leonard examined the local 
enor due to the lad that the whole blob is made to move at the same velocity, whereas in 

leahty the fluid region that coincides with it is strained. He did not emphasize the error 
(iur to tlir initial discretization. 

I lie mathematical papers all show that convergence is obtained when N oo h -> 0 

and - ■ 0, but it not as fast as h (i.e., rr/h • oc ). This last condition was unexpected 

, Leona rd 1980); it means that the vortices must become closely spaced, but also overlap 

more and more for the method to converge. Empirical tests by Nakamura et at. (1982) and 

erlman ( 19, Sb) lully confirm Un- need for the increased overlap. Hald and Beale k Majda 

suggest simple laws like it ('It' 1 with (' some constant and the power i/ less than 1 The 

Mud.es by Hald, Leonard and Beale k Majda all show that the error is proportional to 

som. power ol it, not h. For instance with a Caussian core function f the error is of order 
rr~ . 

I he initialization oi the vortices is more involved than one might expect, even in the 
Ml, , pie ea,e of an unbounded domain. Suppose one has an initial vorticity field u> 0 . One 
lays a grid ol cells oi size h over the domain. One procedure is to place a vortex at the 
ol ,,i:rh °‘ 11, iU,<l 1u !t ,l "‘ circulation within that cell (the integral of u>„ over the 
‘. r '' VU! ' ,HI, - V " Ulroid would then be preferable to the geometric center of the cell. 

Am.Me , ,s to g.ve ,t the value of at the center, multiplied by //-*. The first procedure 
" l, ‘" Ullm,lv< ' ; "" 1 ls preferable if is not smooth; the second one is formally more 
accurate and is recommended if one is seeking better than second-order accuracy and u>„ 
is very smooth. A third procedure is to require that the vorticitv given by (23) equal 

, n,ll< ‘ r 1,1 iht ' K ri<1 l ,uill,s ,,r ««*t of points with a least -squares algorithm 

Nakamura «/ ol. 1982). I his results in a matrix equation for the F/s. This idea is 

attractive, but the implementation can be awkward, and the matrix is very ill-conditioned 
i! rr/h is lur^r. 


Coniu njt tier in practice 

( onvergence has been tested and observed only for unbounded flows (wall-bounded 
Hows are discussed in tj ft.ti). Many papers contain some useful information about the 
behavior ol the method (Nakamura ,1 ol. 19.S2, lb-ale Majda 198.9). but Perlman's 
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1985 study is the most complete. She chose some simple, smooth, radially symmetric 
Euler solutions and systematically computed the errors in actual calculations, including 
the various components introduced by Beale k Majda (1982). She found that for smooth 
u> 0 the method did converge, with the expected order of accuracy. For less-smooth u>„ 
the higher-order cores (ones which satisfy more moment conditions) did not provide any 
advantage. Values of the power q close to 1 (little overlap of the cores) made the calculation 
lose its accuracy in a shorter time. The disorganization of the vortices, compared with 
their regular initial arrangement, was largely responsible for the loss of accuracy. 

The question of disorganization is very relevant., because in practice the vortices always 
look disorganized, see Fig. 1. In all flows of interest, the product of the time and the 
velocity gradients (strain rate, vorticity) is large, so that the fluid gets highly distorted. 
The advantage of the vortex method over the few grid-based Lagrangian methods that have 
been tried is p r ecise!y that it can tolerate this distortion. However, its order of accuracy 
seems to drop to lower values. Quite piobabiy, most vortex calculations produce results 
that are quantitatively accurate (especially the global quantities, e.g., forces and moments) 
but are far out of the range where the theoretical error estimates apply (§/.,?, 9.6). This 
is at the same time a weakness, and a strength. 
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4. Three-dimensional flows" 

I { .i. Essential dijjeriices with two-dimensional flows* 

There are two differences: tin* constraint- on the vortioity to be divergence- free is not 
trivial, and infinitely thin vortices are not a viable approximation. The constraints were 
discussed in jj I./ r In three dimensions an arbitrary collection of Dirac functions as in (17), 
or of their smoothed counterparts as in (23), does not satisfy V.u; 0 and therefore cannot 
represent a vort icily field. This difficulty has been addressed i n two ways. One is to chain 
the vortex elements together to form filaments or nets, which automatically imposes the 
constraint. The other is to ignore the constraint at the level of the discretization and to 
rely on the accuracy of the method: as the method converges the constraint is better and 
better satisfied, as is the other equation (4). The first approach produced the filament and 
segment methods, while the second produced the ‘‘monopole” methods. 

The second difference is also a major difficulty. The velocity of a thin vortex filament of 
circulation 1\ radius of curvature /?, core radius a, and binormal vector b is approximately 
given by 


U 


r 

47 tH 



(30) 


see Hatrhelor 19(17, p. -723. Tims filaments with curvature have a self-induction, and 
moreover it is not well behaved in the limit of an infinitely thin filament (a / li * 0). The 
next term in the approximation, filter (30), strongly depends on the shape of the vorticity 
distribution. <\g.. a ‘‘lop hat” or a (bvussian as in viscous flows (21). This behavior of 
the sell -induced velocity lias nio'ivated a number of studies based or, the Local Induction 
Approximation (LIA) (Hama 1002). (mushier a smooth vortex with circulation V and 
global length scale A, and assume that a/ L -♦ 0. If oik* integrates from time 0 to T and 
t he product /’ V log( /, / a ) / L z is kept constant the motion of t-l.e vortex tends to a finite 
limit as o I. ■ (1, given by (3D). Observe t hat ( 30) depends only on the local characteristics 
of t he vortex. 

I he LIA is reminiscent of the 21) point-vortex approximation since it describes the 
motion in the limit of very thin vortices, Out is in fact almost the opposite since the 
internal dynamics, which were negligible in 21), now control everything. Also, two distinct 
filaments do not interact at all, and different parts of the same filament interact weakly 
through the derivatives that enter the definition of the tangent, the binomial, and the 
radius of curvature. The LIA became especially popular when Hasimoto (1972) showed 
that it could support solitons. However in practice it is not easy to justify the value of a or 
to decide whether a should vary in space or time. The restriction a << L is also unrealistic 
in any situation with viscous or turbulent stresses. At this point this approximation has 
given excellent results in a few specific cases (Leonard but it is not regarded as a 

numerical method adapted to the solution of practical problems. 
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4.2. Filament methods * 

Filament methods have been introduced and refined by Leonard and his co-workers over 
the years (Leonard 1985). Since vorticity is naturally arranged on closed lines Leonard 
chooses a discretization based on vortex filaments each with a constant circulation The 
filaments then induce a velocity field on each other and on themselves, and their motion is 
computed. The self-induced velocity requires special care because the Biot-Savart kernel 
(3) is strongly singular and many different regularizations have been tried. The method 
is in a state of flux, because the regularization that has been used up to now is accurate 
for long waves (it agrees with (30)) but is not satisfactory for short waves on the filament, 
often resulting in a spurious instability (Leonard 1985). Recent work by Winckelmans & 
Leonard (1987) focuses on a solution to this problem which actually borrows a term from 
the LIA and adds it to the usual regularized contribution. 

Leonard’s method is by far the most used in 3D (Leonard 1980, 1985, Ashurst 1983, 
Ashurst & Meiburg 1985, Nakamura, Leonard & Spalart 1986). See also the closely- 
related 3D vortex-in-cell method of Couet, Buneman, & Leonard (1981). With proper 
care this method can be applied to a variety of vorticity-dominated flows. Some involve 
isolated filaments, in which case the core parameters have a st rong influence on the solution; 
others involve “bundles” of numerical filaments to describe a physical vortex, which makes 
the calculation much more expensive but less sensitive to the parameters. A systematic 
application to flows past solid bodies, with boundary layers, separation, and so on, still 
has to he made. 


4.3. Segment methods * 

Chorin (1980, 1982) proposed a method in which chains of vortex elements can also be 
identified, hut they are treated totally separately in terms of the velocity they induce, and 
the simplest and least computationally expensive regularization of the kernel is used (the 
velocity is corrected within a sphere centered on the mid-point of the segment). Chorin 
does not require a single chain to behave like a physical vortex, and relies on bundles or 
clouds of elements to obtain the correct behavior (e.g., (30)). Compared with filament 
methods this method is somewhat cheaper computationally ami offers easier bookkeeping. 
In fact, one could generalize the pattern of elements from a bundle of lines to a “net”. One 
would keep track of a set of nodes, which get moved in tune following the fluid. There 
would also he a set of segments, each with a circulation, a starting node and an end node. 
Any number of segments could start from, or end on, a single node provided that their 
circulations (taken with the appropriate sign) add up to zero. With such a method it 
would he much easier to merge nodes to control their number if needed (§0.^). 

Chorin (1982) used his method to compute the evolution of a disturbed vortex in a 
periodic box, with emphasis on the question of singularities in the 3D Kuler equations 
and the fractal character of the support of vorticity in turbulent flow. He also studied 
boundary-layer transition (Chorin 1980). 

4-4- Monopolt methods* 

A number of three-dimensional i. thuds have been proposed in which the computational 
elements are totally disconnected instead of being arranged in filaments or m*ts (Rehbach 
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1977, Saffman 1980, Beale k Majda 1982, Novikov 1983, C. Greengard 1984, Cantaloube 
k Hubcrson 1984, Anderson k C. Greengard 1985, Mosher 1985). The elements are 
.sometimes called “vortex monopoles” or “vortons”. In such methods the rotation and 
stretching of vorticity by the u>.VU term is explicitly computed instead of being a by- 
product of the relat ive motion of the points that describe the filament. The methods differ 
primarily in the way this Imn is computed. 

I hr attractive features of such methods are: simplicity and computational convenience, 
possibly easier mat hematical analysis, and the freedom for the elements to rearrange them- 
selves so that the topology of the vortex lines changes (the “reconnection process”). The 
objection was raised early on (by the author and certainly by many others) that in such a 
method the vorticity field carried by the computational elements may not be divergence- 
free, which is unnatural. 

Suppose one has a set of elements which carry vorticity in three-dimensional space. Call 
the vector field formed by adding the vorticity of all the elements, as in (17) or (23) but 
with vector values. 1 lie subscript “a 1 stands for “apparent”, as will become clear. The 
procedure is to apply the Biot-Savart law (3) to u; a to obtain the velocity field U. Let us 
define the true vorticity field u;< by — V x U. Clearly the divergence of is 0, since 
it is the curl of some vector field. If the divergence of the apparent vorticity is nonzero, 
:1 < annot be equal to the true vorticity: u) a u? t . More precisely, u > t is the projection of 
u’ (i onto the space of divergence-free vector fields (this is so because the Biot-Savart law 
actually solves the Poisson equation V*U V x u> Q , so that the only link between u? a 
and is that they have the same curl). 

Ci aphically, what this means is that when the “vortex lines” of u; fl terminate within 
the domain (i.e., V\u.> 4 , / 0) the vortex lines of u) t correct tins by connecting to another 
element, or the opposite end of the same element, in as smooth a fashion as possible. 

I bus the small region in which w’ rl / 0 is surrounded by a much larger “halo” in which 
( lb 1 his is disturbing, because one applies the vorticity transport equation (including 
the vorl ex- sf retching and rotation term, etc.) to u; a , not to . Thus if there is rotation 
locally where u/,, is located, the whole halo of rotates too, which ir just not correct. 

If the visualization of a monopole calculation shows that the elements are scrambled 
and some of the vortex lines of u j a terminrte within the domain, then u> a and u? t must be 
strongly different, the calculation :s unreliable, and it should be terminated. Monopole 
methods arc unlikely to be "robust” in the sense of allowing safe integration to large times. 
Kehbach ( 1977) reports serious problems far downstream of bis wing and had to artificially 
reduce the strength oi the vortex elements in that region to obtain a stable solution. 

Nalliuau k w Meiron s ( l9X(i) study is quite critical of vorton methods. They showed that, 
in general, a system of vortons is not a weak solution of (4) and does not even conserve 
tlm integrated vorticity <»r momentum, given by (9) and (10) (see also Leonard 198f>); 
ll»e\ also comment on tlx* central role ol the nonzero divergence of u>„. Nonconservation of 
integrated vorticity is not too serious, because it is the integral of u> a which is not conserved 
always integrates to 0, assuming fast enough decay at large distances). On the other 
hand, nonconservation of momentum is a serious flaw (§7.^). 

1 he most popular application of monopole methods has been to the reconnection pro- 
cess (<\ (Jreengard 1981, Mosher 19Ku), which filament methods are unable to model 
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without outside intervention. Rehbach (1977) also studied the flow past sharp-edged 
three-dimensional wings, and Cantaloube & Huberson (1984) studied propellers and ro- 
tors. Only Greengard had a mathematical convergence proof for his method and attempted 
(with some degree of success) to demonstrate the convergence in actual calculations. The 
three-dimensional problem is probably the most open and challenging in the field of vortex 
methods. 
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5. Extension to viscous and compressible flows 

5.1. Viscous flows 

The original Navier-Stokes equations express local interactions between neighboring 
points, through the spatial derivatives. In contrast, vortices interact at any distance 
through the Biot-Savart law and do not recognize their neighbors. The unfortunate con- 
sequence is that while we were able to formulate one set of equations, the Euler equations, 
very elegantly the necessary spatial derivatives are not available to add any other terms, 
for instance the viscous terms, in a straightforward manner. One needs to alter (19), (20), 
or (23). 

One obvious idea is to use the exact viscous solution (21) to make vortex blobs, at least 
in 2D. The blobs then spread in time like \fvi. Equations (19) and (20) are retained. This 
method is attractive, but inconsistent. In a vortex-blob method, to obtain convergence 
of the transport process one needs to let the core size a tend to 0 (see §3.2). If a is 
made proportional to y/vt , the user can’t control it any more, and the core-deformation 
error cannot be made to vanish. Thus the method “mimics” viscosity, but damages the 
inviscid process. C. Greengard (1985) actually proved that the method solves “the wrong 
equation”: the vorticity is “correctly diffused, but incorrectly convected, even in the limit 
of infinitely many vortices”. 

The most used treatment of viscosity is the random-walk algorithm (Chorin 1973). 
Chorin adds a random component of appropriate magnitude (proportional to \fvKt) to 
the position of each vortex at each time step. In the average, with a large number of 
vortices, this introduces the desired diffusion. When solid bodies are present there is also 
a creation of new vortices along the wall at each step to account for the flux of vorticity 
out of the wall. This approach is very simple, consistent, and partial proofs of convergence 
are starting to appear (Bald 198-4, Goodman 1987). On the other hand it seems inefficient 
to carefully track a vortex if it contains only such a small amount of information (since 
it takes many vortices for the laws of statistics to turn the random walk process into a 
meaningful macroscopic diffusion). Milinazzo & Saffman (1977) created a controversy by 
testing th e method and showing very slow convergence to a simple exact solution; the error 
was of the order of N J /-, with ;V the number of vortices. There are some questions about 
which measure of the error is most meaningful with the random vortex method, but the 
consensus seems to be that slow convergence is observed notably near boundaries (Chorin 
1980), and the ()(N 1/J ) estimate was confirmed by Roberts (1985). 

Probably the most convincing application with random walk was the study by Ghoniem 
k Gagnon (1987) They succeeded in simulating the flow over a backward-facing step 
in the laminar range, up to Reynolds number 229, using about. 1000 vortices. In this 
flow viscous stresses control the flow even away from the wall, so the agreement with 
experiments and other calculations is very encouraging. On the other hand, in the author’s 
estimation results of the same quality could be obtained with a grid-based method with 
much less computational expense, and the competitiveness of the random-walk method in 
terms of computational speed is not excellent . This led the present author to choose a less 
ambitious approach, wherein the viscous effects are totally neglected away from the wall, 
and are included only in a relatively crude manner in the boundary layers. This approach is 
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clearly not universal, but with careful implementation and high enough Reynolds numbers 
it lias been justified by the results obtained on a number of cases. 

Two other possibilities for a viscous method should be mentioned. One could be called 
the “exploding vortex” method: at each step each vortex would be split into a few vortices 
with a spread again proportional to yfcKt. This method leads to a rapid proliferation 
o vortices, and therefore is extremely expensive and impractical. Finally, there is the 
idea of transferring circulation between blobs at each time step in an appropriate manner 
o represent the diffusion (G. Winckelmans, personal communication 1987, after work 
by Cottet and Haviart’s group in France). Thus, (19) is modified. There is no core 
spreading, so the convergence of the invkcid part of the dynamics is not jeopardized. 
C irculation is conserved, but the vorticity centroid apparently is not (one could modify (20) 
to compensate). The method needs to have some “dormant” vortices with zero circulation 
in any region that may become vortical at later times, and the vortices must be closely 
packed, only a distance of order s/uKt apart. In grid-based methods the spacing between 

points is typically of the order of y/cf (where T is a global time scale), which is much 
larger than vvAt. 

The profusion of candidates for the viscous treatment (we even omitted a few that 
showed very little promise) reflects the intense need for such an addition, as well as the 
lack of a truly satisfactory solution to this date. 

5.2. Compressible flows * 

Most flows of interest are at least slightly compressible, and the capability to treat such 
flows with a vortex method would be very valuable. Furthermore one of the key ingredients 
of the method, Kelvin’s theorem of conservation of circulation, still applies to barotropic 
compressible flows. Shocks and other sources of vorticity invalidate the theorem; this 
creation of vorticity could in principle be computed, but this has not been attempted yet. 
bet ns restrict the discussion to shock-free flows. Then ( 19 ) and (20) apply. 

The equation that does not apply any more is (18). The velocity now also depends on 
the dilatation rate («, }- v y ), which is nonzero. One possible approach is to include this 
quantity in the description, and have a number of sources in the flow just like there are 
vortices. T he obstacle here is that the evolution equation for tl. dilatation rate is far from 
being as simple as that for vorticity. Furthermore the dilatation rate, unlike the vorticity 
is not normally confined to a small region. Thus one of the main advantages o' the vortex 
method is lost. These two obstacles seem to have prevented any sustained attempt using 
this approach. The author is only aware of some work along these lines by Professor J. C 
Wu, using a finite-difference method (personal communication, 1983). 

Another approach is to consider flows at low Mach numbers and to approximate the 
dilatation rate. One would then obtain an expansion in terms of the Mach number, similar 
to the well-known and successful Prandtl-Glauert approximation. This was attempted 
by Deflenbaugh & Shivananda (1980). They obtained a Poisson equation for the first 
correction, of the order the Mach number squared, and used a finite-difference method 
U) so vo it. ll IS not dear to the author that the procedure was consistent, because time 
derivatives were found both on the left- and right-hand side of the equations. Indeed the 
numerical solutions were unstable in time, and Deffenbaugh K- Shivananda had to delete 
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some of the time-dependent terms to obtain bounded solutions. They conjectured that a 
more accurate scheme would alleviate the problem. To our knowledge, this approach has 
not been followed upon. 

For future work, it seems that the objections to the first approach are deeper than those 
to the second one. One should be able to formulate a correction at low Mach number. 
However it is not clear that this can be done without the use oi spatial derivatives as in a 
Poisson equation; these can be handled on a grid, but losing the grid-free character of the 
method is a heavy price to pay. The best by far would be to find a grid-free strategy for 
the compressibility correction. 
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6.1. Use of boundary elements 

The traditional way of implementing the inviscid boundary condition (i.e., zero normal 
velocity at a surface) is by using image vortices. However this is possible only when the 
shape of the surface is very simple (e.g., straight line or circle) or when one has an analytical 
conformal mapping from the desired shape to one of the simple shapes. This has led to 
numerous studies of the flow past circles, ellipses, and Joukovsky airfoils, for instance. 
Another serious drawback of images in practice is that a vortex close to the wall is also 
close to its own image. A strong interaction can develop, with the vortex and the image 
forming a pair atid traveling along the wall at a high velocity and in the direction opposite 
to that of the boundary layer. 

In practice the method of boundary elements, in which vortices, sources, or doublets of 
appropriate strength are placed along the boundary to ensure that the boundary condition 
is satisfied, is much more flexible. In other contexts it is called the panel method. At the 
expense of solving a linear system, one can treat, arbitrary shapes. The matrix of the linear 
system is unfortunately full, but in a time-dependent computation it needs to be inverted 
only once so that the cost at each step is of order M 2 and not M 3 (M being the number 
of boundary elements). 

When the boundary-element method is used as part of a vortex algorithm, as opposed 
to a potential-flow solver, it is convenient to use vortices for the boundary elements (Lewis 
1981). Furthermore in some cases using vortices is much more natural (for instance the 
splitter plate for a. mixing layer carries a. velocity jump and therefore circulation). Let us 
restrict our attention to two dimensions for now. Let the points 1 < J < M, describe the 
boundary. The bound vortices are placed at these points. One could construct an algorithm 
to cancel the normal velocity at or near the points xj (for instance at (x, +x J + 1 )/2). 
However, there are large locai velocities and depending on the exact choice of the test 
point and other parameters there is a danger of having fluid “leak” through the boundary 
in the average. It is much preferable to use the stream function: the condition is then 
0(Xjii) - 0(Xj)* Thus we are imposing zero mass flux between Xj and Xj+ j. This 

method is more robust, and also more convenient when one has several surfaces (e.g., the 
walls ol a wind-tunnel ) and wishes to impose the mass flux between them: the condition 
becomes — 0(Xj) Q if x 3 is on one w’all, x ; ^ i is on the other and Q is the 

mass flux. I he stream-function method results in a matrix which is well structured; the 
largest elements are on the main diagonal and the one just below it, and Gauss elimination 
without pivoting is sufficient to solve the system. 

lypically the stream function results from the linear combination of three components. 
If there is flow at infinity this corresponds to a stream function ij' oo(x) - (x x 
where is the uniform freestream velocity. Note that, a slight generalization is possible 
in 21): if the freestream flow has uniform vorticity the vortex method can still be applied 
to the deviation from the freestream vorticity (i.e., Eq. (7) still holds, but (4) does not 
in dl)). I hen there is t he stream function associated with the free vortices, \j'j y which is 
given by (28). The initialization or generation of these vortices strongly depends on the 
problem at hand. Finally there is the stream function of the bound vortices, 0*,, which is 
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given by (28' too. Thus the equations are 

V’6(Xj+l) - V’b(Xj) = V’oo (Xj) - V , oo(Xj + 1 ) + V»/(*j) - ;(xi+i) - (31) 

The unknowns are the circulations F' ; of the M bound vortices. After they have been 
computed, the velocity of the free vortices is computed and they are advanced for one time 
step according to (20). The displacement of the free vortices (and possibly the creation of 
new ones) alters </’/, so that (31) has to be solved again for the next step. 

The question arises immediately that (31) can be written only for j < M; there are 
only M 1 equations, for the M unknowns T^. This is the well-known paradox: the total 
circulation is undetermined. The solution is to add an M th equation that prescribes the 
sum of the circulations. Thus all the elements of the last row of the matrix are 1. Typically, 
one forces the total circulation in the plane to be zero, so that the M ih equation is: 

M N 

En = -£ r ' < 32 > 

i »=i 

(recall that N is the number of free vortices). However this condition is not unique; 
sometimes symmetry or other considerations are applied instead. See the discussion of 
(38) in § 7.2. 

For the extension to three dimensions, one will not have a stream function. One can of 
course revert to a condition on the normal velocity itself, but imposing the condition in an 
integral sense as in (?? ) is preferable. Suppose one has cells on the surface (e.g., triangles); 
the border of each cell defines a contour. One could use a vector potential and require 
that its circulation around each contour is zero. This is equivalent to requiring that the 
(lux of the velocity field through each cell is zero. Since mass conservation is built into 
the discretization, it seems consistent, to use for the boundary an algorithm that exactly 
conserves mass. 

An application of the vortex method with inviscid boundaries is given in §0. Other 
applications are nozzles or wind tunnels. S pal art. (1984) assessed tunnel blockage effects 
l>v imposing the inviscid condition on the tunnel walls (on which separation did not. occur) 
and viscous conditions on the body. One can also represent- a. lifting suiface as a.n inviscid 
boundary (Katz 198|), but then the circulation in the lower and upper boundary layers 
cannot be distinguished, so that one cannot predict viscous effects such as separation; one 
also needs to impose the Kut.ta condition. 

(i. 2. Hi mark on the choice of two numerical parameters * 

in general, there are no explicit equations linking the different length parameters, so 
that the user is left with his/her intuition to choose values. However in the simple rase 
of an inviscid boundary one can use the following simple argument to relate the spacing 
between the points Xj and x M i and the core radius <r. We idealize the boundary as a 
straight line on the .r axis with equally distributed points x, at intervals A. bet V denote 
the velocity jump across the sheet. The ideal stream function is t />, V\y\/2\ the flow is 

uniform at ) V f'l above and below the sheet. With the discretization by vortices, the flow 
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is nearly uniform away from the sheet by forms the usual “cat’s eyes” around the vortices. 
The stream function discretized using vortex blobs with the core function F == r 2 /(r 2 + 1) 
(see (46)), is given by 


V’rf(x) - - 


UA 

4tt 


E lo «( 


j~ oo 


Ix-jAe^ + ^x 
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(33) 


1 he circulation of each vortex is — U A, and the denominator inside the logarithm represents 
an additive constant that was set to ensure that at x - Xj = j Ae x , V’i = = 0 (recall that 

the stream function will be sampled at the points Xj to write (31)). One can verify that 
for large \y\ the stream function r(> d given by (33) becomes arbitrarily close to -U\y\/2 + C> 
where C is a constant. The mass flux between the test point Xj and a point far up will 
be correct only if C — 0. This mass flux is an important quantity when the walls are the 
walls of a tunnel or nozzle, especially of variable section. It is less important in external 
flows. A numerical summation of the series (33) showed that (7 = 0 is obtained when < 7 / A 
is equal to 0.153. The optimum value of <j/A naturally depends on the core profile used 
(the F function). In practice, the value 0.153 can rarely be exactly respected since the 
spacing of the points is generally uneven, but it is a useful guide. 


6.3. Example: a curved mixing layer 

Mixing layers, either temporally- or spatially-developing, have been a frequent subject 
for vortex methods (Ashurst 1979, Aref k Siggia 1980, Ashurst & Meiburg 1985). A short, 
unpublished study of curved, spatially-developing mixing layers will be described here 
primarily as an example of how realistic boundary conditions can be applied. The study 
was suggested by Dr. N. Mansour (NASA Ames Research Center) after some experimental 
studies showed a strong effect, of curvature on the structure of mixing layers (with the faster 
stream inside, the layer was much more disorganized and three-dimensional). 

Although this was not the only way to proceed (see the treatment of the flat mixing 
layer using two semi-infinite vortex sheets in Leonard 1980) it was decided to use curved 
walls as guides for the flow. Let R be the radius of the mixing layer, R x the radius of 
the inner wall and R 2 that of the outer wall. The values chosen satisfied R^/R r~ 0.9 and 
R'lj R 1.1, and the “test section” covered one radian; see Fig. 3. Boundary conditions 
were enforced on the two walls, on an inflow boundary and on a section of the splitter plate 
using bound vortices. On the walls and plate, the inviscid no-penetration condition was 
applied (most- mixing-layer simulations do not enforce any condition on the plate; they just 
place a row of vortices of equal circulation on it). Along the inflow boundary, the velocity 
normal to the boundary («#) was imposed and equal to the velocity of the irrotational flow 
(he., proportional to 1 /r). This is an example of (31) with Qj ^ 0. It can be interpreted 
as a porous boundary with imposed blowing. The system was closed by requiring that 
the circulation of the bound and the free vortices add up to zero, see (32). This value 
was arbitrary but ensured that the flow outside the walls was relatively quiescent, so that 
the velocity field was smooth at the inflow and outflow boundaries. Provided that the 
test section is long enough, the manner in which the linear system is closed (the M ih 
equation ) is unimportant, as are the inflow and outflow conditions; here we are invoking 
Saint Ven a nl 's principle. 
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Figure 3. Computation of curved mixing layers, a) Fast fluid outside, U 2 ~ 2U\\ 
l>) fast fluid inside, V 2 -■ V l / 2. o free vortices; + vortices bound to 
walls; A vortices bound to splitter plate; □ vortices bound to inflow 
boundary 
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No condition was needed at the outflow. The trailing edge of the plate gives an example 
of a smooth transition from the hound vortices to the free vortices; the two families of 
vortices have the same spacing, core size, and essentially the same circulation (this assumes 
the proper relationship between the mean velocity, the time step, and the spacing). At 
every time step the last bound vortex becomes a free vortex, keeping its circulation. The 
vortex that is farthest downstream is deleted. Thus the calculation can be run as long as 
desired. 

In parallel with the numerical study, the linear stability of curved vortex sheets was 
investigated. By generalizing the well-known derivation for flat sheets (Batchelor 1967, p. 
511) it was found that the growth rate a for a wave with wavenumber a is given by: 

» 2 = j(", -(M 2 + - up (34) 

where V x is the velocity of the inside stream and U 2 is the velocity of the outside stream. 
Note that this formula reduces to Batchelor’s result (a — a\U x U 2 1/2) as aR — » oo, that 
is, as the radius of curvature becomes much larger than the wavelength. 

Equation (54) shows that the mixing layer is more unstable if the inner stream is faster 
{I J \ > However the relative difference in growth rates is smallest for the most 

unstable waves (largest |a|) so this effect is not significant. Figure 3 shows mixing layers 
with velocity ratios V 2 /V x of 2 and 1/2 respectively, and confirms that even a significant 
curvature (compared with the size of the largest eddies) has little effect on two-dimensional 
mixing layers. This strongly suggests that the third dimension is necessary to explain the 
effect observed in experiments. 

Another point of interest is that the boundary condition on the splitter plate did not 
seem to make a noticeable difference in the flow pattern. There was a possibility that the 
pressure disturbances from downstream would influence the circulation of the new vortex 
that is created at each step, thus creating a feedback mechanism with the sharp trailing 
edge as an amplifier. However, in the calculations this eflecl was very weak, and the flow 
did not differ noticeably when the boundary condition was enforced using Ml) and when 
the vortices hound to the plate wen* given inform circulations. 
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7. Viscous boundaries 

7.1. Application of the Biot-Savart law inside a solid body 

In §6 we showed how the inviscid, no- penetration condition can be conveniently enforced 
using boundary elements, preferably vortices. We are now interested in satisfying also the 
no-slip condition. A priori this is more difficult, but we have a very fortunate result that 
shows that, for external flow's, this can still be achieved with a single set of boundary 
elements (Wu & Sankar 1980). Thus, calculations with both the no-penetration and no- 
slip conditions arc almost as simple as those with only the former. The key idea is to 
extend the velocity field and the stream function into the region occupied by the body. 

Consider a solid body occupying the bounded region D , with boundary dD. It may be 
multiply-connected, made up of sub-regions Dj. For now we restrict ourselves to a fixed 
body and two dimensions. Suppose we obtained a vorticity field such that the stream 
function it induces through the Biot-Savart law (8) satisfies the no-penetration condition, 
dll' /'Os — 0 on dD. Recall that (8) represents the Green’s function for the whole plane. 
We may have achieved this using boundary elements; in any case there are only vorticity 
elements (so that i /> is defined and single- valued), and none of them are inside D. The 
formula (8) can be applied for x inside D. and there V 2 V> = u> — 0. Therefore V’ satisfies: 

^ 0 in D, (35a) 

0 cm dD . (356) 

as 

This is Laplace's equation with a Diricblet boumiary condition. It is a well-posed problem 
with a unique solution up to an additive constant in each subregion: t/’ ~ 

Therefore U, as Riven by (fi), is identically zero in D , Also consider cty/fJn, the normal 
derivative of ty as one approaches ()!? from the inside. Since ip is constant in D ls 

-r- ■ 0 on dD, (3G) 

On 

i.e., the tangential velocity just below the boundary is zero. This is a very useful result. 

Several remarks. First, one could have first imposed d^/dn -- 0 and obtained dif'/ds - 0 
as the by-product (Laplace's equation with Neumann condition also has a unique solution 
up to an additive constant). In other words, if (35a) holds (35b) and (3(>) are equivalent. In 
the numerical method we elected to explicitly apply the no-penetration condition d^/ds - 
0 because it recpiires less smoothness of 0 (observe that (31) is written directly in terms 
of </\ not in terms of derivatives) and also because it seemed more important to conserve 
mass exactly; the viscous effects, including the no- slip condition, are treated a little more 
loosely. Second, the solution to (35) is unique only if D is hounded and dD is closed, 
that is, for an external flow. Third, the argument on the normal derivative applies only if 
one can approach dD from the inside, that is, if D has a finite thickness. The argument 
totally breaks down for an infinitely thin body, e.g.. a p’ate, and in th.? numerical context 
one needs to be careft.l if the t hickness of the body is not large enough compared with the 
spacing of the vortices along the wall. 
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Figure 4a. Flow pa.s1 a bluff body, just after impulsive start. — boundary; 
• vortices; — — velocity vectors 


The behavior of the velocity field is illustrated in Fig. 4. A simple body shape was 
considered, just after an impulsive start (i.e., there are no free vortices yet). Thus we have 
irrotational flow around the body. Boundary vortices were used to enforce the difi/ds = 0 
condition, discretized as in (31). Tlu* velocity vectors computed using (26) are plotted 
at the nodes of a regular grid, including inside the body. In Fig. 4a one recognizes the 
usual potential-flow solution outside the body. Ins»de the body, the velocity vectors are 
very small, showing that the iormal argument base i on Laplace’s equation (35) does work 
well nun erically. This was not entirely obvious, because usually the function 6„ in (23) 
does not have bounded support, so that there is a little residual vorticity inside the body, 
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and also because we forced t/’ to take the same value (C\) only at the discr te points x ; 
(see (31)). The natural smoothness of solutions to Laplace’s equation helps. In Fig. 4 the 
boundary elements represent a vortex sheet, and the strength dT/ds of the sheet gives the 
velocity jump across the sheet. Therefore dT/ds gives the tangential velocity above the 
elements. Thus once the boundary vortices are known it is a very simple matter to obtain 
the limiting velocity of the potential flow at the wall, and the surface pressure through 
Rernoulli’s equation. Thus even for potential flow the method of boundary vortices has a 
distinct advantage over those ba... d wn sources or doublets. 

7.2, Creation of vortieity 

The interpretation now differs for inviscid and viscous flow. In viscous flow the vortices 
along the boundary are not considered as mathematical tools to enforce a boundary con- 
dition, but as actual vortieity that is imparted to the fluid by the shear stress at the wall. 
Solid walls have long been interpreted as sources of vortieity (Lighthill 1963). New vortices 
are created at every step and released into the fluid, just like at the trailing edge for the 
mixing layer (§tf.3), but now all along the wall (Oborin 1973). The circulation of the new 
vortices is computed as is they were bound vortices, but they are immediately considered 
as free vortices and allowed to move with the fluid. Symbolically, their contribution to ip 
is transferred from \ to </'/, on the right-hand side of (31 ). so that t/’oo + V’/ alone satisfies 
the boundary condition. The vortices then move according to (20), by an amount that is 
()(At). This alters -I iJ'j by an amount 0(A/), so when (31) is solved again tpb and 
therefore the circulation of the new vortices are 0( At). Thus, a small amount of vortieity 
is released at each step. The vortices are created not right on the boundary, but a small 
distance </ 0 outside it. There are “wall points” x ; where (31) (with Qj - 0) is applied, 
and “creation points”, say x' where the new vortices are placed. 

The wall points and creation points can be distinguished in Fig. 4b, which is a detail 
of Fig. 4a. Notice the smooth potential flow away from the wall, the velocity vectors 
becoming very small inside the body, ami the transition region in which the direction of 
t he velocity fluctuates. The microscopic description of the wall (on the scale of \Xj+ \ - x ; | ) 
is a little crude, but on a macroscopic level the method has the following advantage. The 
row of vortices in Fig. 4b (and in general with any attached boundary layer) forms a vortex 
sheet . The tangential velocity above the shear layer is f^ t ., the velocity of ihe potential 
flow’, and the velocity below it is zero. Therefore the vortices translate along *he wall at 
a velocity V r j 2. This is known to be correct, in the sense that in a boundary layer the 
average velocity of the vortieity, defined by 



is precisely equal to l\./2 (just use w ~u fc and an integration by parts). The flux of 
vortieity along the wall is an important quantity and is correct (as mentioned earlier, it 
can be very inaccurate when image vortices are used). This is another example of how 
small scale errors can cancel each ot her so that integrated quantities ( here, U av ) are correct. 

As in the inviscid Hows, we nerd an extra equation to close the system (31) for earh 
body. This condition is that the circulations of the vortices created along each body <t 
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Figure 4b, Flow past a bluff body, just after impulsive start (detail of Fig. 4a). 

-o- boundary and wail points x^; - vortices at creation points x^; 

i 


velocity vectors 


each time step add up to zero (Wu k. Sankar 1980): 




Note that tins condition is rigorous, unlike (32) which looks the same but was arbitrary 
(actually, (38) was known and (32) was chosen by analogy). If one considers the whole flow, 
(38) implies that the total circulation is conserved, a well-known result, see (13), However, 
(38) is a stronger condition since it applies to each body separately. In summary, for each 
body described by M wall points there are M 1 type-(31 ) equations and one of type (38). 
Note that the type-(3l) equations are coupled between bodies, so that the whole matrix 
is full (except for the few type-(38) lines). 

In inviscid flow, w can have a singularity of the order of a one-dimensional Dirac dis- 
tribution, and V’ is continuous but with a jump of its normal derivative across i)D . In 
viscous flow, has a jump across 01) but is hounded, and both derivatives of t/’ are con- 
tinuous. Therefore (3G) applies on either side of 01),, i.e,, the no-slip condition is satisfied. 
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u> is bounded because with viscosity it cannot remain in a sheet; instead it diffuses into 
the fluid. The mechanism which prevents it from accumulating at the wall is viscous; 
thus the vortex-creation algorithm is consistent only with a viscous flow solver. In actual 
calculations numerical effects play a strong role in this process, see §7.^ and S.2. 

Several slight generalizations are possible. If the body is translating, this 'an be taken 
into account using the Q } ' s in (31). If the body is rotating at an angular velocity Q the 
situation is a little more involved since there is vorticity u> •- 2fl inside D. Again, one 
adjusts the Q } \. The stream function induced by the inside vorticity can be computed 
with the Contour Dynamics formula (Zabusky, Hughes, k Roberts 1979) and subtracted 
from the right-hand side of (31). This was done by Spalart k Leonard (1981) for the 
dynamic-stall study, see Fig. 7. On the other hand, if the body is deforming or several 
bodies are in relative motion there is an extra difficulty because the matrix involved in 
(31) becomes time-dependent, so that it has to be computed and inverted at each step. 
Finally, in three dimensions the same argument as in (35) can be applied to ‘he vector 
potential. Thus the extension to 31) is possible, but it has not been made yet; we first 
need to obtain a robust method for unbounded 3D flows. 


7 . ('hoivt of numerical parameters m 

In tin* creation algorithm just dcs<rihed, three length scales are involved near the wall. 
As before (tj fi.2) we have A, the spacing between the wall points, and <r, the core radius. 
Now there is also </ 0 , the distance from the wall to the creation points. The requirements 
conflict. d ( , should be large enough comparer! wiih o for the residual vorticitv inside 
the body to be small, but it should also be small enough compared with A to maintain 
a strong coupling between the creation point and the wall point (so that the matrix is 
well-conditioned). Moreover, a should l>c large enough compared with A for the cores to 
overlap in the s direction, so that the velocity field is as smooth as possible (fig. 4b). The 
cliflioi.lt y arises because the vortices are circular in a region, the boundary layer, where the 
length scales in the .« and n direction are naturally disproportionate. This is partly what 
motivates separate treatments of the boundary layer (see j t 8). 

l or viscous flows, a quantitative mnditioi as in has not been obtained; t/v cannot be 
forced to give the right mass flux from Xj; to oo without making the vortex cores impinge 
on the boundary, in that sense, the actual boundary is blurred, somewhere between the 
wall points and the creation points. Figure 4b shows a satisfactory set of parameters: 
d 0 A,, /?, tr - A u /4, where A„ is an average over the boundary: A - S/M where S is 
the length of the boundary. The need for a balance between A, c/ 0 and o implies that the 
local value of A should not have very wide variations. It is helpful to cluster the points 
(thus reducing A) in the regions with high curvature and near sharp corners, bur this 
clustering should be moderate. Typically, A should not vary by more than a factor of 2. 


Kutta condition 

The Kutta condition, like conformal mappings and imsge vortices, is a major tool of 
classical hydrodynamic theory, und it has been explicitly incorporated in many vortex 
methods (e.g.. Stansby ,V Dixon (19X2)). Actually, a generalization of the classical condi 
lion (which just provides the circulation for steady flow on an airfoil with a single sharp 
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trailing edge) was needed since the flows have vorticity and are time-dependent, and there 
may be several sharp edges, like on a square. However the idea is the same: the velocity 
at any sharp corner must be finite. Note that “Kutta conditions” have been applied to 
smooth surfaces at separation, which is probably an improper use of the wcia (in fact, one 
is just applying the no-slip condition). Note also that a complete viscous method does not 
need the Tutta condition, because iliis condition is simply the result of separation of the 
boundary layers, which an accurate viscous calculation can account for. 

It turns out that the method outlined in §7.2, with creation of vorticity and the no-slip 
condition satisfied all along dl), does a good job of reproducing this viscous behavior, even 
when the viscous terms are not included. For instance in Fig. 1 the flow near the various 
trailing edges and convex sharp corners is physically correct, although no special treatment 
or extra condition was applied there. When confronted with a convex sharp corner, the 
method naturally expels boundary-layer vorticity into the outer how until a smooth flow 
without large velocities around the corner is established. The high curvature causes a rapid 
acceleration and deceleration of the boundary layer, which causes separation. Vorticity is 
released until the streamline originating at the corner is aligned with the upstream side of 
the corner. 

There are two major reasons for this “pseudo-viscous” behavior. The first is that the 
vortices are created at a small but nonzero distance do from the wall so that whenever 
there is a deceleration of the boundary layer they acquire a small velocity component 
directed away from the wall, which is what initiates separation. The second is that time- 
integration errors, especially with large velocities ?md strong streamline curvature, move 
the vortices away from the wall. The analogies between integration errors and viscous 
effects are discussed in more detail in §<S\2 and 9.L In any case this somewhat fortuitous 
behavior is consistent, and the flow near sharp edges has always been physically correct. 
Note that the shedding of vorticity can be only transient as on a starling airfoil without 
stall, or permanent as in Fig. 1 or in the case of flow past a square. The separation from 
smooth surfaces is a much more delicate effect, and is discussed in detail in §8. 

7.7. Example: starling vortex on an airfoil 

The flow development around an airfoil at incidence started impulsively from rest is 
considered, and will help clarify the issue of the Kutta condition. This is a classical 
problem (Prandtl k Tietjcns 1934). The airfoil is NACA 0012; the integral boundary- 
layer treatment is activated (see §<S\i?) and the Reynolds number is high. Figure 5 shows 
four snapshots of the flow. At first one sees a sheet of vortices leaving the leading edge and 
curling up. The starting vortex then grows in size. Since total circulation is conserved, an 
opposite circulation remains on the airfoil, contained in the boundary layers. A sheet of 
vort ices connects the primary vortex and the trailing edge. This sheet clearly has a nonzero 
circulation (velocity jump), since it undergoes a Kelvin Helmholtz instability. Thus, the 
circulation on the airfoil is progressively brought up to its final value. In fact, the downwash 
of the vortex on the airfoil reduces the effective angle of attack, and tends to 0 only like 
1 ft. The last snapshot is at much larger time; the wake now has zero circulation (note 
that it does not roll up any more) and the flow leaves the trailing edge smoothly. The 
force vector at early times shows a lower lift than in the final state, as well as some drag 
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Figure 5 Airfoil starting impulsively without initial circulation. 

a) iUoo/c = 0.02; b) tU^/c = 0.2; c) tU^'c = 0.3; d) W^/c - oo. 
c is the chord and the velocity. • • ■ vortices; I force 


(induced drag from the starting vortex). At large times the drag is very small. 

7.6. Pressure and force extraction * 

The pressure was eliminated when the curl of (lb) was ♦ dcen to derive (4) and is not 
necessary to advance the vorticity equation, but one often needs it as a diagnostic or for 
the treatment of boundary layers. The wall pressure and the forces on the bodies can be 
obtained directly from the rate of creation of vorticity. Consider a wall with s and n the 
tangential and normal coordinates. Equation (7) can be written u, = -V.(Uu> - i/Vu>), 
i.e., the flux of vorticity is Uw - iA7u>, which reduces to -uVu at the wall since U = 0. We 
take the dot product with the normal vector n to obtain the flux through the boundary: 
-vdv/dn. This is the rate at which the wall is “creating vorticity”. Now consider the 
momentum equation (lb); at the wall it reduces to Vp = ^V 2 U or equivalently Vp = 
x u>) (using a vector identity). In 2D, or in 3D with a flat surface, one can show 
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that n x (V x u>) + n.(Vw) - 0 (in 3D, curvature seems to introduce another term). In 
those cases we then get n x Vp - vdw/dn, or in 2D 

dp du> __ o 2 r 

ds V dn ~ dsdt ' 

The notation d 2 Y/dsdt is appropriate since this is a rate of creation of circulation per unit 
length and unit time. This quantity is known, since we are creating the new vortices at 
every time step. It is then a simple matter to obtain the wall pressure in 2D: (39) gives us 

the derivative of the pressure along the wall. We usually apply the trapezoidal rule and 
write 

^ (40) 

where the new vortices of circulation T' and F' + , are created during a step of length At. 
If we integrate (39) around the wall, since p is single- valued we get 

/ d 2 Y 

}ala! d ^ 0 ' <«> 

that is, the total circulation created along each body is 0. This is how (38) was derived 
by Wu & Sankar (1980). I hus if we apply (38) the wall pressure given by (39) returns 
to exactly the same value after going around the body, which is not always the case with 
other methods for computing the pressure. A weakness of (39) is that it gives p only up 
to an additive constant. In many cases, the front part of the body is in contact with 
irrotational fluid from upstream; one can then obtain a good estimate of the additive 
constant by assigning Die stagnation pressure to the front stagnation point (this neglects 
the unsteadiness in Bernoulli’s equation). Also, Ibis method does not provide the pressure 
away from the wall. One could however use this wall pressure as boundary condition 

for the Poisson equation V 2 p = - V.(U.VU) (although the additive constants would be 
trou blosome). 

The equation dp/d. s - d 2 l'/0sdt which led to (40) is much less sensitive to the viscous 
details of the boundary layer (more specifically, the exact distribution of the vorticity in the 
normal direction) than dp/ds ----- uduj/di, (used by Chorin (1973)). Thus (40) is very well 
adapted to vortex methods and an accurate wall pressure can be obtained even with a crude 
representation of the boundary layer itself. For instance if one has a steady attached flow 
the flux of vorticity along the wall is It; ,2 as mentioned earlier (§7.2). It is balanced by a 
flux of vorticity through the wall at a rate d(U 2 /2)/ds; thus we get dp/ds = -d(U 2 /2)/ds, 
i.e., the well-known Bernoulli equation. Many procedures have been proposed to obtain 
the pressure from the time-dependent version of Bernoulli’s equation; these methods are 
very cumbersome in separated flows because of the numerous branch cuts in the velocity 
potential (one for each vortex). 

The global pressure force on a body is given by 


e. x 


(42) 
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After integration by parts and use of (39) this becomes 

e * * / EFt x ds ’ (43) 

which is the negative of (14), applied to the newly-created vorticity. This is not a coin- 
cidence. Neglect the viscosity for now. What happens is that the motion of the vortices 
under (26) preserves (14) (the momentum of the fluid), and that the creation of vorticity 
effects the exchange of momentum between body and fluid. Thus (43) amounts to a global 
momentum balance (but it cl be applied to each body separately). One can check that 
(43) gives the right “apparent mass” for a body undergoing acceleration. A similar formula 
holds for the moment on the body, using the invariant (15) (Wu & Sankar 1980). If the 
viscosity is not neglected, for instance if a random walk is used, then (14) is not (and 
should not be) conserved when a vortex collides with the wall, and this is how the friction 
force comes in. In practice, even without random walk there are such collisions due to 
noise in the velocity field near the boundary (see § 8.2 and 9.3) so that the pressure force 
and the momentum balance do not exactly agree, but this is a small effect, and (39) and 
(43) have been very satisfactory at high Reynolds numbers. 
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8. Separate treatment of the boundary layers 

8.1. Overview 

The boundary layers are the most difficult regions for a vortex method for at least three 
reasons. First, the standard method uses radially symmetric elements and therefore is not 
adapted to the disparate scales near the wall (with finite-difference methods, high-aspect- 
ratio grids are routinely used near the wall). This results ui significant noise in the normal 
component of velocity; recall Fig. 4b. Second, near the wall the viscous terms become 
dominant, and the available vortex methods do not treat them very efficiently. Third, the 
boundary-layer region is usually close to being steady in Eulerian coordinates: du/dt is 
small, in contrast with the outer region where it is the Langrangian derivative Du/ Dt 
which is small. Thus Eulerian methods have an intrinsic advantage in the boundary layer. 

One response to the problem of disparate scales is Teng’s (1982) method, which uses 
elliptic vortices. Teng applied it only to a boundary layer, using identical ellipses and a 
random walk for the viscosity. He obtained good results, but apparently the method has 
not been used further. One could certainly extend the idea by using ellipses of different 
eccentricities; the vortices would be very flat near the wall and circular away from it. The 
ellipse eccentricity and orientation would be simple functions of the distance to the nearest 
wall. Note that such a concept is quite different from one in which the ellipse parameters 
would obey dynamical equations reflecting the straining and the self-induced rotation. The 
ellipses are introduced for kinematic reasons only; to obtain a more realistic flow field near 
the wall. A drawback of Teng’s method is that the formula for the velocity induced by 
an elliptic vortex is much more complex than for circular blobs, so that the method can 
be competitive only if an ellipse can replace at teast several blobs. It would pay to find a 
simpler formula for an “elliptical blob”. 

Many methods can be found in the literature in which the no-slip condition is not 
enforced using the boundary elements as described in §7. Most often the no-penetration 
condition is enforced using images, and there are a few “separation points ’ at which vortices 
are created and released into the flow (Lewis 198.1, Stansby A' Dixon 1982, etc.). When 
the body has convex sharp corners, they are obviously going to cause separation, and it 
is consistent to apply a Kutta condition. The rate of release of vorticity is ffr/2, with 
the appropriate sign, where (' , tne edge velocity of the upstream boundary layer. This 
prescribes the circulation of the new vortex, ar.c. the Kutta condition provides a condition 
on its position. Such methods satisfy the intuition and have the major advantage that 
only a few new vortices are created at every step, compared with hundreds in the method 
presented in §7.2. They are typical of applications on minicomputers. On the other hand, 
the use of the Kutta condition is still very delicate. Another weakness is that separation, 
in the sense of a transfer of vorticity from the wall region to the outer region, does not 
occur only at t he obvious points. “Secondary separation” occurs at random on the back 
face of bluff bodies, and its neglect may be largely responsible for the excess circulation 
found by some studies and sometimes compensated for by arbitrary means (see § 9 . 6 ) . 
Stansby & Dixon (1982) emphasize the importance of secondary separation. In the long 
run, and especially for smooth bodies, methods based on a few selected separation points 
will probably be superseded by those w.ncu allow separation in any region, depending on 
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the boundary-layer dynamics. 

The idea of using an inviscid solver for the outer flow and a boundary-layer solver near 
the wall is old, and has been very successful for attached flows. In such a case one is 
primarily interested in the viscous friction drag, and possibly in the displacement effect, 
on the outer flow. The trouble starts when the flow separates; the standard boundary- 
ayer equations develop the CJoldstein singularity. Any proposal to use the boundary-layer 
equations for a separated flow should include a discussion of this singularity. It is all too 
easy to mistake its eilects for purely numerical problems, or conversely to suppress them 
using numerical dissipation. 


S-t. Simple delay of separation 


As mentioned already the vortices induce a fair amount of noise near the wall, which 
tends to scatter them. As a result, some move away from the wall, and others collide with 
it. 1 he ones that collide can be pushed back, or they can be absorbed and their circulation 
be reintroduced a distance d (l from the wall (see §9.3). Statistically, the vorticity tends to 
f n t away from the wall. As soon as there is an adverse pressure gradient, the straining 
adds to this effect, and separation occurs. This is a mild case of the situation described for 
sharp edges in §7./,. The pure vortex method consistently predicts separation at the onset, 
of an adverse pressure gradient, which is essentially the behavior of a laminar boundary 
layer. With turbulent boundary layers, separation can take place much farther along the 
wall and a 21) method is unable to reproduce this by itself. 

This led to the idea of predicting the separation point by a separate method, which can 
ta < t ic Reynolds number and the effects of i urbulencc- into account, and intervening in 
t ic dynamics of the vortices to delay their separation up to the correct position. If this 
inti 1 \en1 ion is p< 1 foi nu <1 in a sensible manner, we can retain the smooth blending between 
the at tuchcd boundary-layer vorticity and the outside vorticity (this is not achieved with 

the methods that release vortices only at a few points and do not enforce the no-slip 
condition). 


At every time step, the pressure distribution is computed using (40). With it., the two 
boundary layers starting at the front stagnation point of each body are calc ulated. So far 
only integral methods fThwaites' and Head's, see Oebeci K' Bradshaw 1977) have been used,’ 
m conjunction with CJranville's transition criterion, and the time-dependent effects on the 
boundary layers have been neglected. Thus, there is room for experimentation with more 
sophisticated boundary-layer solvers. The only information needed from the boundary- 
layer solution are the separation points (we do not. attempt, to continue the boundary-layer 
solution beyond separation). On each body, downstream of these two primary separation 

points the vortices are five- to behave normally, so that secondary separation can take 
place. 

Upstream of the separation points indicated by the boundary layer solution, we prevent 
the statistical drift, and separation of the vortices simply by re-absorbing the new vortices 
after only one step. A fresli regular row of vortices carries the boundary- layer vorticity as in 
hg. 4, representing a thm attached boundary layer. Just beyond the indicated separation 
point, the vortices are suddenly free and subjected to the straining associated with the 
adverse pressure gradient, and they leave the vicinity of the wall. This procedure achieves 
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Figure 6. Simulation of the flow over an airfoil, a) without delay of separation; 
b) with delay. • • • vortices; f lift and drag 


the goal of making the vortices separate approximately in the right place. 

Figure 6 shows the flow over an airfoil calculated without and with separation delay. 
In Fig. Ca, observe how the vortices remain attached on the lower surface, thanks to the 
favorable pressure gradient, but separate just downstream of the leading edge on the upper 
surface, in a mild adverse gradient. This is what laminar boundary layers would do. In 
Fig. 6b with a very high Reynolds number the integral boundary-layer solvers predicted 
turbulence, and separation only at the trailing edge. This activated the separation-delay 
device and we are left with ait ached boundary layers and a thin, smooth wake. The lift is 
higher, and there is almost no drag. Separation delay was also applied in Fig. 1; it made 
a difference on the upper surface ol the slat, and on element 4. It was fully responsible for 
the Reynolds- number dependence that was described (§/./). 

This device has been used extensively, and is very helpful in many cases (Fig. 1, 5, 6b, 
7, 13 and 14). It. has weaknesses. It cannot account for reattaclmient. The weakest part 
of the solution is the transition prediction (but this is true of any method). To obtain 
a flow pattern that agreed with experiment in Fig. 7, the transition criterion had to be 
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adjusted; the increase in R e from the instability point to the transition point was taken 
as half of Granville’s correlation. Another minor problem occurs when one cannot find 
the front stagnation point because the flow is too disturbed in front of the body. This 
occurs in the rotating-stall .al dilation, see Fig. 14. In that case, one just bypasses the 
boundary-layer solution for this body at this step. Finally, the integral methods have a 
tendency to predict separation a few wall points downstream of sharp corners, and one 
needs to override this prediction and set the vortices free a few points upstream of the 
corner. Because oi these issues, when starting with a new shape it is good practice to first 
simulate the flow without separation delay, then activate the device if there is a need for 
it, and carefully check the behavior of the boundary layers. 

Figure 7 shows a case in which separation delay made a significant difference. An NACA 
0012 airfoil is oscillating from 5° to 25“ at a reduced frequency of 0.25 (Spalart k Leonard 
1981 ). The overall flow strongly depends upon separation on the upper surface just behind 
the leading edge, and this separation is controlled by transition in the boundary layer. The 
Reynolds number was 2.5 10 6 . The forces and moments on the airfoil were in very good 
the agreement with experiments (Spalart, Leonard k Baganoff 1983). 


8.3. Chorm’s tile method* 

In 1978 Chorin introduced a vortex-like method for the boundary layer. As before, the 
elements carry circulation, move with the fluid, and undergo a random walk for viscous 
diffusion. However they are now "sheets” or “tiles” and the velocity is computed in a much 
simpler manner thanks to the boundary-layer approximations (the tiles interact only with 
their neighbors). This method is much cheaper computationally than the usual method 
or l<ng s method (1982), and is designed for use in a hybrid algorithm, with conventional 
vorticts away fiom the wall. ( ireulation is exchanged by the two subdomains, primarily 
in the direction sheet > blob. Chorin discusses ihe Goldstein singularity and the need 
foi a viscotis-inviscid coupling strategy to remove it, but the subsequent users of the tile 
method completely overlooked this issue. 

C heer ( 1983) applied the complete hybrid method to the flow past a circle and Joukovsky 
airfoils, with Reynolds numbers of about 1000. She obtained flows with “laminar” behavior, 
i.e., separation as soon as the boundary layer enters the adverse pressure gradient, and 
satisfactory drag and lilt values. Ghoniem k Gagnon (19*7) applied it to an internal flow 
at rather low Reynolds numbers, which seems a little superfluous since there were no thin 
boundary layers (the inflow was a parabolic profile); quite probably a pure vortex-blob 
method would do just as well. In summary, the tile method succeeded in introducing 
the boundary-layer character for the near wall flow; it is fast and solves the first problem 
described at the beginning of t f 8 (the scale disparity), but not really the other two. The 
convergence of the random- walk method is rather slow here as it is with blobs. Also, the 
Goldstein singularity is a major problem and Chorin himself (1978) did not assert that it 
was fully solved. 


8./,. Coupling with a finite -difference solver* 

As mentioned above, buleriun methods have significant advantages for boundary layers, 
especially if they are steady or quasi-steady, which makes the time integration undemand- 
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ing. Furthermore, the location of the steep gradients is known and it is easy to cluster 
the grid points near the wall. This motivated an effort to couple an inviscid vortex solver, 
for the outer region, and a finile-oiLerence solver, for the wall region (Spalart, Leonard & 
BaganofF 1983). 



Figure 8. Sketch of the hybrid finite-dilference/vortex algorithm, o vortices; — 
grid lilies; — • velocity 


The inner region is a thin strip along the wall, see Fig. 8. The boundary-layer vorticity 
equations are solved using a standard, implicit, centered linite-ditlerence scheme (T. Pul- 
liam, personal communication 1981 ). Simple transition and turbulence models are used. 
The challenge is in the coupling with the outer region. The two regions exchange circu- 
lation as in the tile method, according to the direction of the velocity at the interface. 
The coupling of the velocity fields is more delicate. A viscous-inviscid coupling was imple- 
mented based on the displacemc 1 thickness (the centroid of the inner-region vorticity). 
The procedure is intuitively correct and very similar to the ones used in attached flows, 
but there is no proof that it actually removes the Goldstein singularity, especially when 
large amounts of vorticity c ross the boundary. Some numerical difficulties (oscillations) in 
the separation region suggest that it does not (see Spalart et al. 1983 for details). Partly 
for this reason, this method has not been as robust, and has received much less use than 
the simpler “delay" method 

The method however gave encouraging results for the flow past a circle. The “drag 
crisis’ 1 in this flow at Reynolds numbers near 3 • 10 s is one of the major challenges in 
the field of viscous flows. The flow was calculated foi Reynolds numbers between 10 4 
(the lower limb for thin boundary layers to form) to 3 > 10* (beyond which there were 
no experiments). Figure 9 shows that a decrease in drag coeflicient was observed, but 
that the crisis itself (including ihe sudden drop to very low drag values and the steady 
lifting states observed in experiments) was not. The drag den eased at higher Reynolds 
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Figure 9, Drag coefficient of a circular cylinder. — average of experiments; 
o calculation 


numbers because the boundary layers transitioned; the turbulence model then kept the 
vorticity attached longer, so that it crossed the boundary into the outer region farther 
downstream. This is illustrated in Fig. 10; the mean streamlines, shown at Re — 10 5 and 
10 6 , reflect the position of the separation region. By symmetry, only half of the flow needs 
to be shown. The Reynolds-number effect was qualitatively correct, and the quantitative 
agreement with experiments, over a wide range of Reynolds numbers, was fair. 

For the future, the use of an Eulerian solver for the inner region should be the most 
efficient procedure. On the other hand, it is not the most elegant since it is a hybrid, 
and the coupling and singularity issues definitely need to be re-examined. One could of 
course solve the full Navier-Stokes equations in the inn'T region; the obstacle there is 
that one would nerd to solve a Poisson equation to obtain the velocity field, instead of just 
integrating the equations = -u> and v y =- -u r up from the wall. The velocity boundary 
condition on the interface would be a serious problem. 
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9. Practical aspects 


9.1. Time -integration schemes 

A wide variety of time-integration schemes ace used in Computational Fluid Dynamics, 
and the choice is often a matter of taste in addition to the technical corcerns. Still, one 
can narrow down the choice considerably in the case of vortex methods. The first point is 
that all the schemes used have been explicit schemes. To use an implicit scheme one would 
need to invert (26) or at least a linearized form of it, and (26) is too complex for this to be 
done efficiently. Moreover the usual stability advantage of implicit methods is much less 
crucial with vortex methods than with grid-based methods (see below). The second point 
is that the most costly part of the algorithm, by far, is the evaluation of the time derivative 
using (26). Also, the calculations use a relatively small amount of memory, so that storing 
old values of the time derivative is not a problem. This suggests that multi-step schemes 
are preferable to predictor-corrector schemes. 

Based on these remarks the author has used the Adams- Bash forth second-order scheme, 


x,(/ i A/)* Xj{t) 


A/ V2 U(Xj, 0~ 2 U(x,,/ 


AO 


(44) 


for the integration of (20), which represents a good compromise between accuracy, sim- 
plicity ami cost. As usual, the first step for each new vortex is done using the Euler 
scheme. A majority of the studies in the literature used either the Euler first order scheme 
or Kunge Kutta schemes of second or fourth order (the Runge-Kutta second-order scheme 
is also known as “modified Euler” or “Huen’s scheme” ). The use of the low-order Euler 
scheme, when the Adams- Bashforth scheme is hardly more costly, was motivated either by 
a question of consistency with the random- walk algorithm (in (Turin's method) or by a 
faulty concept of numerical dissipation, as discussed below. 

The integration of (20) is most, difficult when large accelerations are experienced by a 
vortex, typically when it is close to another vortex. Therefore a good model problem is 
the motion of just two vortices (Spalart k Leonard 19*1). The problem can be expressed 
in complex notation and simplified to 


dz til 

<1 It 

:(*>) 1 


(45a) 
(45 h) 


where ; is the separation vector between the two vortices, :* is its conjugate and il is a 
frequency. The exact solution is : «xp(iQ/): the vortices orbit around each other with 

angular velocity il. Their distance is constant. It is a good exercise to program (45) and 
observe the behavior with different numerical schemes. 

One finds that with any scheme, if the time step A t is large compared with 1 /il the 
solution is inaccurate. However due to its nonlinearity (45a) reacts very differently from 
the linear equations that are usually studied, e.g., the equation dz/dt ill: whiih lias 
the same exact solution, i'sually, as an explicit scheme becomes inaccurate (because A / 
is too large) it also becomes unstable and the solution grows to infinity (it “blows up'). 
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With (45a) the d' lance between the vortices grows until \z\~/ti is of the order of A/, 
then settles down. The vo. tires now orbit at a lower angular velocity, Q/\z\ 2 . With some 
schemes \z\ will shrink again and in some cases one gets an oscillation, especially with 
multistep schemes as the primary and spurious roots alternatively dominate (Spalart & 
Leonard 1981). With the Adams- B ash forth scheme, such unphysical behavior has not been 
observed. 

The lesson to be learned here, which holds in practice with many vortices and has 
been recognized by many authors (Milinazzo fc Saffinan 1977, etc.), is that in almost 
all situations time-integration errors tend to scatter the vortices. The scattering occurs 
whenever the local rotation rate 12 is large compared witli 1/A/ and corrects itself, since 
the scattering precisely reduces L2. The corollary is that a sustained instability, with the 
solution “blowing up” to infinity, never occurs. In that sense, the method is “robust". 
However, a solution can remain bounded and be very inaccurate. Also, a scattering of 
vorticity reduces the velocity fluctuations and the kinetic energy. In that sense, it represents 
a numerical dissipation. Lagrangian methods, unlike fixed-grid Eulerian methods, allow a 
perfect convection of flow structures by a uniform velocity field; however when gradients 
are present inaccuracies show up as with any other method. 

The scattering of vortices has implications for the invariants (11, 12, 15, 16). The total 
vorticity (9, 13) and the impulse (10, 14) are linear quantities in terms of the vortex 
positions Xj. Therefore with the common time-integration schemes which are all linear, 
if (9, 10, 13, 14) are conserved by the spatial discretization they are still conserved even 
allowing for time- oitegrat ion errors in (20). This is not so for the nonlinear invariants 
(11, 12, 15, 16); they are only "semi-conserved 1 ' and have often been used as measures of 
the integration errors. For evample Nakamura, Leonard N' Spalart ( 1982) defined a global 
measure of the numerical dissipation, an "effective viscosity”, »ased the time evolution 
of (15), prin inrily to verify i hat it was purely caused by t iwi< -integ ation errors (indeed, 
it was). Some researchers pushed tins reasoning one step too far, concluding that time- 
integration errors were desirable, since they introduce the viscosity that is so d.flvult 
to introduce by other means! i hey then use the most inaccurate and dissipative (when 
applied to (45)) scheme, the first-order Kuhn scheme. Nakamura it al. never implied that 
the elective viscosity had any meaning locally. Moreover, if errors are what one is after, 
it is cheaper to just increase the time step, rather than switching to a lower-order scheme. 

One final remark about steady flows. If the velocity field is naturally time-independent, 
grid based methods can usually find the solution faster by direct solution or by relaxation 
feclm. ties than by a straightforward integration in time. Vortex methods do not have 
this option and treat any flow as time-dependent; the vortices move even in a steady How. 
fortunately, most sepaiated flows arc strongly unsteady, so this difference is not of much 
importance in practice. 


.V.J. f‘im functions 

In a sepurated-How calculation one will need to compute /’ and F many times, so that 
it is desirable for both to have simple formulas. ( ‘oust quent ly it is convenient to choose F 
first. As mentioned earlier the (iaussiau core (21; (but with a fixed radius <r substituted 
for v'W) i s attractive hut is expensive to compute, and yields only second-order accuracy. 
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A very simple core profile was used in all the cases shown in these notes: 


F(r) = 


r 

7-2 + 1 ' 


( 46 ) 


]1 is second-order accurate too. Its only weakness is that the vorticity decays rather 
slowly at large distances, like r~ 4 , so that the residual vorticity mentioned in §7 .3. is not 
negligible. A vorticity distribution with bounded support or at least faster than r 4 decay 
would be preferable as far as the interaction with walls is concerned. Hald (1979) and 
Beale &; Majda (1982, 1985) give formulas for higher-order cores. Recall that Perlman 
(1985) found that these performed better only if the vortices remained well-ordered, so 
that in complex separated flows their advantage is probably marginal. 


9.3. Control of the vortex count 

In most applications, many new vortices are created at every time step and added to 
the list. For instance in the calculation of Fig. 1, there are about 1300 vortices in total, 
and 249 new ones are created along the wall at each step. Clearly, useful calculations are 
possible only if vortices are deleted at about, the same rate they are created. 

First take the simple case of the mixing layer of Fig. 2. At each step, one new vortex 
is created at the trailing edge, and one “old” vortex is deleted near the outflow boundary. 
The easiest procedure is to delete the oldest vortex; a slightly better one is to delete the 
vort ex that is fart hest downst ream. This makes little difference, because the last 10 or 20% 
of the domain are not included in the “measurements” (e.g., velocities, Reynolds stresses) 
anyway since that region is disturbed by the lack of vortices beyond the outflow boundary. 
Either method allows one to generate a satisfactory statistically-steady state of the flow, 
and to maintain it as long as desired. 

In the less simple cases with solid bodies, two operations allow one to control the number 
of vortices. The fiist. one is .he absorption of vortices by the wall. When a vortex collides 
with the wall, one has several options. One can reflect it as in an elastic collision, or 
just push it back to the wall. A third option which is consistent with the treatment of 
viscous boundaries described in §7.2 is simply to delete the vortex. If this is done befoie 
the evaluation of </>/ for (31 ), the right-hand side of (31) will reflect the loss of circulation 
near the wall, and the new vortices will make up for it. (also adjust the last line of the 
right-hand side, (32)). This absorption procedure has the advantage of helping to contain 
the number of vortices, and of slightly reducing the noise ne .r the wall. 

The absorption of vortices is helpful, but is not enough since it is not under the user’s 
control; a majoiity of vortices still leave the vicinity of the wall and move into the wake. 
Thus one needs a more active and adjustable means of limiting N. This can be done by 
merging vortices. The method used is adapted from Rogallo’s (personal communication, 
NASA Ames Research Center. 1979). Deffenbaugh k Marshall (197(i) also used a simple 
merging device. At every step, pairs of vortices are merged into single vortices. The goal 
is to achieve this with the least amount of modification of the velocity field. 

Since a large local modification is unavoidable, it makes sense to require that the velocity 
in the far field be altered as litth- as possible. Suppose the vortices j and k are merged, 
and the new vortex is given the circulation F' and the position Using complex notation 
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and omitting sonic constant factors, tlic change in the velocity at 


The Taylor expansion for large jej is 


+ 0 ( 1 - 


I he two leading terms are removed by taking V 1 - Vj -j V h and z 1 (L — + r***)/r # : 
the new vortex inherits t he sum of the circulations of the old ones, and is placed at their 
centroid. This could he expected; the circulation (13) and the vortex centroid (14) are 
conserved. On the other hand, (15) and (16) cannot he conserved (recall that they were 
already not conserved because of time-integration errors). 

I he next term in (48) cannot he removed and we take it as the merging criterion, i.e., 
we merge a pair only if this term is below some tolerance. After insertion of T' and J and 
taking of the modulus it becomes: 

iiji'.i i~i - 


Consider the two factors in thus est imat e. The first, one encourages merging of weak vortices 
(small r j and/or I 1 *), which is sensible since weak vortices cost as much as strong ones, 
but carry less information. It also discourages merging of vortices with nearly opposite 
circulations (small \Vj I l\-\) which also makes sense since in such a merging z f ends up 
far away from z } and zt? which is not natural. 

1 he second factor encourages merging of ( lose vortices (small \z } c^j) as one would 
expect. We still need to prescribe j;| in the denominator The most sensitive place in the 
flow is the wall region ; therefore we insert the distances an 1 d ^ from the vortices to the 
closest wall. The formula that was finally adopted is 




a,)V-> 


Notice how tli#' formula was symmetrized with respect to dj and </^.. The length parameter 
!) yt :uiovv one to shift resolution closer to the walls (small I)[ j) or into the wake (larger /)(>), 
hut the calculations shove little sensitivity to over a wide range (Spalart, Leonard & 
Ha.ga.nofl 1983). A typical value is 10% of the chord. 

finally, the quantity l 0 in (50) has the dimensions of a velocity and is the tolerance. 
I ypica.1 values are small: less than 10 4 times the frecstrcam velocity. In practice if is most 
convenient to let t he program adjust \ \ f during the run. One chooses a target number of 
vortices, and employs a feedback mechanism that raises V (t if there are too many vortices, 
and vice-versa. The merging procedure described here has been very convenient in practice, 
is intuitively correct and lias a rational basis in tin' Taylor expansions. On the other hand, 
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merging is a very noisy operation locally, and its effect on the convergence is unknown. 
Therefore, it is not appropriate for all applications. 


9-4- Efficient programming: operation counts 

Vortex codes are relatively easy to optimize, because almost all the operations are spent 
computing the N 2 interactions in (26). Thus one really needs to optimize only about ten 
lines of code. There are a number of nearly-obvious tips. The division by 2ir and the x 
product should be done outside the 0{N J ) loop. The formulas for the two interactions, 
k — ► J and j — > k, have a lot in common, namely tne computation of 


(Xj - x k ) 


P 




x*l\ 


(51) 


After this, the circulations 1\. and i } are presumably different and the overlap ends. One 
should have an outer loop (k - 2,N) and an inner loop ( j = ],fc - 1) and compute (51) 
only once. There is no need to store it (which would occupy N 2 words of memory); use it 
right away to increment U* and U } (or store ail the j k terms and sum them up outside 
the j loop but inside the k loop). In 2D there is no need to compute any square root, 
and one can avoid exponentials. The operation count is 8N 2 with the simplest core, (46), 
and would be hardly higher for more complex but algebraic cores. In 3D, the symmetry 
between k —■ j and j -- k is useful too. One should also precompute quantities like the 
tangent vector outside the 0(N 2 ) loop. 

In practice, periodic calculations could be much more expensive than unbounded ones, 
because the stream function and its derivatives involve the computation of several com- 
plex exponentials N~ times (see §10.1), and exponentials cost many floating-point op- 
erations. However, this drawback can be completely eliminated with a little extra pro- 
gramming. One needs the values of exp(*7r(c* - :j)/p) for all values of k and j. What 
one does before entering the 0{N 2 ) loop of the program is to precompute and store 
exp(*7r ;*.//>) and exp (~ inz h /p) for every k. Then the operation that is done N 2 times 
is exp(/7r;*./p) x exp( - inzj/p): the exponentials have been eliminated. This can easily 
speed up the execution by a factor of 4 or more, depending on the computer. In fact, 
11k i ) ( A ) component of the operation count js only 18N* so a well-programmed periodic 
melho.i is only about twice as expensive as the unbounded-flow method. Thus there is no 
reason to approximate (52) by taking the first few terms in the series, as has sometimes 
been done in t he past. I he same trick applies to the computation of the stream function 
1/7 at the wall points, for (31 ). 


Anothei tip (oinerns the use of complex variables, in 2D. They are convenient, but 
degrade the efficiency of some computers to some extent (the machines that slow down 
when the loops have skips ol 2). Some machines even refuse to vectorize any loops involving 
complex variables. 1 his should be checked carefully before writing a code. 


9.0. Flow charts 

Here we give flow charts of two programs, named KPDl and KPD2. Both are available 
from COSMIC (a part of NASA) in the U. S., as are their periodic counterparts for 
cascades. A listing of KPDl also is included in the author’s thesis (Spalart, Leonard k 
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Baganoff 1983) (there is an error on line 15, p. 107; replace NW by NWALL(b))- KPD1 
is the basic bluff-body code with no boundary layer treatment; KPD2 has the integral- 
method treatment of the boundary layers (the “separation delay” of §8.2). The flow charts 
provide an overview of the material covered until now, and specify the order in which the 
various operations are performed. 



Figure 11. Flow chart of the KPD1 code 


9.6. Assessment of the accuracy 

This can be delicate in practice. Some numerical methods (in particular with explicit 
time-integration schemes) give the user a strong warning when the time step is too long: 
they “blow up”. Other methods (spectral, for instance) generate “wiggles” which are also 
a warning that finer resolution is needed. Vortex methods do not give such warnings; in 
most applications the vortices look disordered, but the calculation never blows up. This 
is partly due to its automatic conservation of some key integral quantities (9, 10, 13, 14). 
Again, the possibility that the calculation is not accurate on small scales but is accurate 
on the intermediate and large scales is likely. 

The most prominent result of a bluff-body calculation is the average drag coefficient 
Co- We treated four shapes with sharp corners, so that the primary separation points 
were fixed, and compared with experimental results (Spalart 1984, and unpublished). The 
shapes were a square, an equilateral triangle in two positions and a vertical flat plate or 
thin airfoil. They were all treated with the same code and roughly the same resolution 
(typically 1000 vortices). The results are summarized in a table (the — * indicates the 
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Figure 12. Flow chart of the KPD2 code 
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The uncertainty is at least ±0.2 on all the values. Even then, the disagreements are 
significant. Note also that the calculation often over-predicts Cd, but not always. In 
general, the shedding frequency (Strouhal number) is better predicted than the drag. These 
results are not very satisfactory, but we were unable to improve them or to convincingly 
demonstrate the convergence of the method, even via large variations of the numerical 
parameters (N, At, d 0 , etc.). The above table provides a severe but conservative estimate 
of the performance of the author’s method. In some cases, the agreement with experiment 
is much better, within 5% (see §10.1). 
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Many of the comparable methods also overpredirt the drag. One reaction to this is to 
introduce a “fudge factor” to bring the drag down. Arguing that cancellation of vorticity 
of opposite signs occurs in reality but not in the calculation (which is not even obvious) 
some authors arbitrarily make the circulation of the vortices decay as they move into the 
wake. By adjusting the decay rate one can get any drag one wants. In the author’s mind 
the idea of invoking Kelvin's theorem, building a method around it, and then violating it 
just to fine-tune one global quantity, C'p< is strange. 

Milinazzo & Saffnian ( 1 977 ) already pointed out that “comparison with gross experimen- 
tal features, such as drag, is not conclusive, and in any case the real flows are turbulent, and 
therefore three-dimensional, at large Reynolds numbers.” This point is valid, especially 
if the experimental flow has large-scale three-dimensionality. It is possible that much of 
the drag difference is due to this effect. The author’s calculation of dynamic stall (Fig. 7) 
agreed with experiment, better than those for steady bodies. In this case, the whole span 
of the wing oscillates in phase, which makes the flow more two dimensional. The degree of 
three-dimensionality of the flow past two-dimensional shapes, such as a circular cylinder 
or a backward-facing step, is a matter of debate and the subject of ongoing experiments. 
Reliable numerical simulations should also be obtained in the next few years. 

In summary, there is no recipe for the estimation of the errors in a vortex calculation of 
a complex flow like Fig. 1. It is essential to monitor the vortices and velocity field as in 
Fig. 1 or lb, but this is only as good as the user’s intuition of what the flow should look 
like. Finally, it should be remembered that we are looking at very challenging flows. If one 
tackled the multi-element airfoil with an existing finite-difference method, very probably 
1 ,ie best one could do is a rather crude calculation which would not enter its convergence 
pattern, at least not with a reasonable number of points. One would also be dependent on 
crude transition and turbulence models similar to t .hose we used, so that the viscous effects 
are very delicate with any method. When applied to complex flows, the vortex method 
will not yield perfect results, but it will yield very useful ones. 
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10. Treatment of cascades and screens 


10.1. Spatially -periodic flows 

The computation of two-dimensional flows that are spatially periodic in one direction is 
almost as simple as that of unbounded flows, and is useful for the treatment of temporally- 
developing mixing layers and of cascades, for instance. Complex notation is convenient 
here. Lamb (1945, p. 224) gives the formula for the stream function generated by a row of 
vortices of circulation V at the points z 0 + np, n being any integer ( z 0 and p are complex): 


V’U) - 





i 

r 


(52) 


Compare this formula with (8). This function can be regularized by a “core” as in un- 
bounded flows and a convenient one is the analog of (46): 



The rest of the method is identical to the unbounded case; one simply substitutes (53), for 
the stream function, and its derivatives, for the velocity. See §0.^ for a tip in programming 
(53). To make the visualizations consistent one chooses a “main period” between some 
value z\ and z\ + p; if a vortex moves out of the main period (i.e., if the real part of 
(z - Z\ )/p leaves the interval [0, lj) the vortex is translated by ±p so that all the vortices 
can be seen next to each other. 

An example of a periodic calculation is the study of a row of chevrons, shown in Fig. 13. 
These chevrons were candidates for an acoustic barrier around the inlet of a wind tunnel at 
NASA Ames Research Center. The pressure drop across the barrier was estimated to be of 
the order of 5 to 10 times the dynamic pressure, based on the ratio of the open area to the 
pitch of the cascade (essentially the argument of Batchelor, 1967, p. 375). The calculation 
predicted a much larger value: 86 times the dynamic pressure. The issue was settled 
when experiments by Professor j. Foss (Michigan State University) produced a value of 
82. Such losses were unacceptable, and the design was dropped. The open-area argument 
vastly underestimated the loss of energy in the flow due to the repeated separations. This 
was one of the more successful and directly useful applications of the method. 


10.2. A simple model for screens* 

Screens are often used in wind tunnels to generate homogeneous turbulence, to introduce 
a pressure drop (e.g., for mixing-layer studies), or sometimes to smooth out an irregular 
flow and prevent separation on a cascade. Here we present, a simple, macroscopic model 
that can account for the latter three effects and is a good exercise in itself. It was used in 
the chevron calculations (see Fig. 13). it could be used for nonperiodic geometries too. 

Consider for simplicity a screen along the y axis in two dimensions. We model it as 
a pressure jump Ap(y) ~ p(0 4 ,y) — p(0“,y). The effect of the screen on the fluid is 
equivalent to a body force f -- Ap(/y )<S(.r )e x with b the one-dimensional Dirac distribution 
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Figure 14. Propagation of a stall cell in a cascade. •** vortices; j forces; 
— streamlines 


the points that describe the screen. For large values of the product (in the limit C — * oo 
one actually has a solid wall) one will need a fully coupled algorithm, instead of one which 
goes back and forth between the fluid and the screen. 

10.3. Compulation of rotating stall in a cascade* 

The periodic version of the method was also applied to two-dimensional cascades of 
airfoils with the hope of observing “rotating stall”, or more appropriately “traveling stall” 
(since the cascade is actually not circular). When a compressor stalls, the stall pattern is 
often uneven from blade to blade, ana stall cells tend to propagate around the compressor 
at a fraction of the velocity of the blades. Such a flow involves a multiply-connected domain 
(since several blades are involved) and large regions of separation, which makes it a good 
candidate for the vortex method. 

Figure 14 is taken from Spalart (1984, i98o). Five NACA 0012 airfoils are included 
within the numerical period p. They are at 45° to the plane of the cascade. The incoming 
flow, of magnitude ( -y, gives them a 30° incidence. Let e denote the chord of the airfoils; 
the pitch was chosen equal to c. The time interval between frames in the figure is 2c/ 
The computation was started from irrotational flow with a dipole in the wake (to break the 
blade- to- blade periodicity) and has been running long enough to reach a well-developed 
state (see Spalart 198n for the transient behavior and other values for the stagger and 
inflow angles). One observes a stall ceil which alters the flow near two of the blades, 
and drastically reduces or even reverses the mass flux between them. The cell propagates 
upwards at about 40% of / V This figure, and the flow pattern, are in qualitative agreement 
with flow visualizations, hut more detailec comparisons were not attempted. The method 
is limited by its reliance on rather crude transition and turbulence models (as in the 
dynamic-stall study), by the resolution wlrch is becoming marginal flOO wall points and 
300 vortices per blade), and by the fact that it is two-dimensional. However Spcziale, Sisto 
& Jonnavithula (1986) extended the study to cambered airfoils and a wider parameter 
range. 
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11. Time-saving algorithms’ 

11.1 . VortfT'iri’Cf ll methods * 

These methods discard the Green Vfunction approach (3, 6, 8) and solve tin* Poisson 
equation for the stream function hy a grid method (Christiansen 1973). Typically, with 
N 9 grid points the solution is obtained in (){ N 9 log( A ' ) ) operations, so if N* is comparable 
with N the vortex-in-cell method is much faster than the ()(N 2 ) method. In practice N ' 
is at least as large as N , and could be much lar ger if the vorticity distribution is very 
intermittent. The procedure is to define a grid and to obtain vorticity values at the grid 
points according to the circulation of the vortices (if any) present in the adjoining cells. The 
vorticity is “transferred 1 ' to the grid, hut only for the Poisson solution; the vortices retain 
their identity. The beauty of this method is that it still has the mw- dispersion advantage of 
other vortex methods. Once the stream function is available it is differentiated to obtain 
the grid values of the velocity. The velocity is then interpolated to the position of the 
vortices and the vortices are moved, thus completing the work for this time step. 

The best applications are those in which A 1 can be kept to a value close to N , that is, 
when the vortical region is rat her predictable and of simple shape, for instance an internal 
flow or a mixing layer (A ref & Siggia 1980, Court, Buneman K r Leonard 1981). In such 
cases, one can save a large fraction of the computational effort. The method is not grid- 
free, and since boundary conditions are needed on V’ the grid has to be body-fitted. This 
would be a serious drawback in engineering applications. Note also that the transfer and 
interpolation steps have a cost only ( )(N ), but are hardly vectorizable. Another topic of 
research is in the local errors associated with the exchange of information between the 
regularly-arranged grid points and the vortices, which can fall in any region of the grid 
cell. The grid destroys the isotropy of the algorithm, and m general a vortex will have a 
spurious self induced velocity. See the work by Conet (i al. (1981) and Anderson (1987)1)). 

11.2. l.umpitHj methods* 

These methods still use the Biot-Savart approach, but reduce the operation count hy 
approximating some of the interaction terms in (2b). using laylor expansions. Spalart K r 
Leonard (1981) presented such an algorit hm, claiming an operation count of order 
L. Creengard Rokhlin ( 1987) claim ()( A ) for their algorithm. Divide the domain where 
tie* vortices are into a regular array of ceils, j he cells are of size /, large enough to contain 
many vortices each (in contrast with the vortex-in-cell method). I hey can overlap with 
1 lie solid body; in that sense the method is “practically 11 grid-free since the grid is not 
body-fitted. Now consider two cells that are separated by at least a few times /. The 
interactions are proportional to 1 / z. , using complex notation. I his function is analytic and 
has a rapidly converging laylor expansion. II Zp and Z/v are the cell centers, z t is in cell 
/., and Zj in cell A, we Taylor-expand the function l/' in tin 1 vicinity of Z/ v Zp\ see 

f ig- h r >- 

There are several options. I he requirement- is to incur a relative error of less than some 
small number t on any ol the terms. Spalart Leonard chose to use a varying number of 
terms depending on tin* value of \Zp Zp jl (more terms ii it is small). If the cells are 
neighbors the interactions are computed vortex by vortex, because the laylor expansion 
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Figure 15. Schematic of cell-to-ccll interactions in the “lumping” methods 


fails (i.e., it is too close to its radius oi convergence). One can optimize the number of 
cells, and the best value is of order Vne overall cost, for a fixed value of t, is then 

of order A r3//2 . Greengard Sc Rokhiin (1987) carry the idea further by nesting the cells. 
There arc many levels of lumping, in essence they keep \Z K - Zi\/l roughly constant, 
instead of varying the number of terms. Spalart Sc Leonard wasted terms on the very - well 
separated cells (large | Zjc — Z/J/7). The nesting removes the need for computing many 
vortex-by-vortex terms. On the other hand, one needs many more cells. Greengard Sc 
Rokhlin’s arrive at an operation count of order JVlog*(c). 

Greengard & Rokhlin’s statement that the cost is of order N is not rigorously correct, 
because e is not a constant (there is no reason why it should equal the precision of the 
computer). For the method to converge, t must tend to zero like some power of A^ as N 
tends to oo. Thus the cost is really of order fVlog 2 (iV). Similarly, Spalart & Leonard’s 
N 3 ' 1 estimate should be revised upwards, at leasv to A’ 3 / 2 log( A f ), because as « — * 0 one 
needs to include more Taylor terms for a given separation. A factor log(A r ) is not very 
important especially since larger values of N improve the throughput of the computer, but 
one likes to be rigorous if possible. 

In practice, Spalart & Leonard’s algorithm was helpful on serial computers (CDC 7600) 
with c = 10 2 , but was hardly worth the effort on vector computers (Cray IS), because it 
makes vectorizatiou difficult and results in short loops. One probably needs thousands oi 
vortices for the lumping to bring a significant improveir.. ut (the gather-scatter capabilities 
of the newest machines may make a difference). Greengard Rokhiin report an impressive 
advantage starting with a few hundred vortices and with < ss 10~ 5 , again on a serial 
machine (a VAX 8600). The performance of their method on a vector machine will be 
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extremely interesting. The nested method is very promising, and the fact the grid does 
not interfere with the body is a definite advantage over vortex-in-cell methods. Also, empty 
cells do not consume any work. The extension to 3D will require more work since we are 
not just considering the function I/ 2 , but is possible. 
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